DUDLEY? 

3RARY 

NAVAL  r 

-  UATESCHOOl 

<ONTE 

jo943-5101 

SECURITY  CLASSIFICATION  OF  THIS  PAGE 


REPORT  DOCUMENTATION  PAGE 


la  REPORT  SECURITY  CLASSIFICATION 
Unclassified 


2a  SECURITY  CLASSIFICATION  AUTHORITY 


2b  DECLASSIFICATION/DOWNGRADING  SCHEDULE 


1b  RESTRICTIVE  MARKINGS 


3  DISTRIBUTION/AVAILABILITY  OF  REPORT 
Approved  for  public  release,  distribution  is  unlimited. 


4  PERFORMING  ORGANIZATION  REPORT  NUMBER(S) 


5  MONITORING  ORGANIZATION  REPORT  NUMBER(S) 


6a   NAME  OF  PERFORMING  ORGANIZATION 
Naval  Postgraduate  School 


6b  OFFICE  SYMBOL 
(If  applicable) 
AA 


7a  NAME  OF  MONITORING  ORGANIZATION 
Naval  Postgraduate  School 


6c  ADDRESS  (C/ty,  State,  and ZIPCode) 
Monterey,  CA  93943-5000 


7b  ADDRESS  (City,  State,  and  ZIP  Code) 
Monterey,  CA  93943-5000 


8a  NAME  OF  FUNDING/SPONSORING 
ORGANIZATION 


8b  OFFICE  SYMBOL 
(If  applicable) 


9  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 


8c  ADDRESS  (City,  State,  and  ZIP  Code) 


10  SOURCE  OF  FUNDING  NUMBERS 


Program  Flement  No 


Project  No 


ldlk    fM^ 


Work  unit  A^ession 
Number 


1 1   TITLE  (Include  Security  Classification) 

Implementation  of  a  Personal  Computer  Based  Parameter  Estimation  Program 


12  PERSONAL  AUTHOR(S)   Graham,  Robert  Gordon 


13a  TYPE  OF  REPORT 
Engineer's  Thesis 


13b  TIME  COVERED 
From  To 


1 4  DATE  OF  REPORT  (year,  month,  day) 
1992  March 


15  PAGE  COUNT 
147 


16  SUPPLEMENTARY  NOTATION 

The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not  reflect  the  official  policy  or  position  of  the  Department  of  Defense  or  the  U.S. 

Government. 


17  COSATI  CODES 


FIELD 


GROUP 


SUBGROUP 


1 8  SUBJECT  TERMS  (continue  on  reverse  if  necessary  and  identify  by  block  number) 
Parameter  estimation 


19  ABSTRACT  (continue  on  reverse  if  necessary  and  identify  by  block  number) 

Aircraft  parameter  estimation  is  the  process  of  extracting  numerical  values  for  aerodynamic  stability  and  control  derivatives  from  flight-test 
time  history  data.  This  process  can  be  used  as  a  verification  or  validation  tool  for  results  obtained  from  wind-tunnel  testing  or  through 
computational  analysis,  and  can  obtain  or  improve  estimations  of  dynamic  derivatives.  This  study  implements  the  MAT1*AK  Personal  Computer 
(PC  1  based  maximum  likelihood  estimation  routine  for  aircraft  longitudinal  and  lateral-directional  derivatives.  The  parameter  estimation  was 
first  accomplished  on  generated  simulated  data,  with  and  without  noise.  The  noise  consisted  of  measurement  and  state  noise  which  used  the 
DrydenGust  Model.   Secondly,  two  actual  longitudinal  flight  test  maneuvers  are  analyzed  for  the  F14A  and  the  T-37  aircraft.  Additionally,  the 
simulated  portion  of  this  study  can  be  an  excellent  instructional  aid  in  Flight  Dynamics  and  Flight  Test  Courses. 


20  DISTRIBUTION/AVAILABILITY  OF  ABSTRACT 

J  UNCLASSIFIED/UNLIMITED  J   SAME  AS  KEPOK1 


D1IC  USERS 


22a  NAME  OF  RESPONSIBLE  INDIVIDUAL 
Richard  M.  Howard 


21   ABSTRACT  SECURITY  CLASSIFICATION 
Unclassified 


22b  TELEPHONE  (Include  Area  code) 
1408)646  2870 


22c  OFFICE  SYMBOL 
A  A;  Ho 


DD  FORM  1473,  84  MAR 


83  APR  edition  may  be  used  until  exhausted 
All  other  editions  are  obsolete 


SECURITY  CLASSIFICATION  OF  THIS  PAGE 

Unclassified 


TP6044S 


Approved  for  public  release;  distribution  is  unlimited. 

IMPLEMENTATION 

OF  A 

PERSONAL  COMPUTER  BASED 

PARAMETER  ESTIMATION 

PROGRAM 

by 

Robert  Gordon  Graham 

Lieutenant  Commander,  United  States  Navy 

B.S.,  United  States  Naval  Academy,  1979 

M.S.,  Naval  Postgraduate  School,  1991 

Submitted  in  partial  fulfillment 
of  the  requirements  for  the  degree  of 

AERONAUTICAL  AND  ASTRONAUTICAL  ENGINEER 

from  the 

NAVAL  POSTGRADUATE  SCHOOL 
r  March,  1992 


ABSTRACT 

Aircraft  parameter  estimation  is  the  process  of  extracting 
numerical  values  for  aerodynamic  stability  and  control  derivatives 
from  flight -test  time  history  data.  This  process  can  be  used  as  a 
verification  or  validation  tool  for  results  obtained  from  wind- 
tunnel  testing  or  through  computational  analysis,  and  can  obtain  or 
improve  estimations  of  dynamic  derivatives. 

This  study  implements  the  MATLAB  Personal  Computer  (PC)  based 
maximum  likelihood  estimation  routine  for  aircraft  longitudinal  and 
lateral-directional  derivatives.  The  parameter  estimation  was 
first  accomplished  on  generated  simulated  data,  with  and  without 
noise.  The  noise  consisted  of  measurement  and  state  noise  which 
used  the  Dryden  Gust  Model.  Secondly,  two  actual  longitudinal 
flight-test  maneuvers  are  analyzed  for  the  F-14A  and  the  T-37 
aircraft.  Additionally,  the  simulated  portion  of  this  study  can  be 
an  excellent  instructional  aid  in  Flight  Dynamics  and  Flight  Test 
Courses . 
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I .   INTRODUCTION 

Aircraft  parameter  estimation  is  the  process  of  extracting 
numerical  values  for  aerodynamic  stability  and  control 
derivatives  from  flight -test  time  history  data.  Aircraft 
flight  tests  designed  for  this  purpose  are  generally  motivated 
by  one  or  a  combination  of  the  following  objectives: 


1.  The  desire  to  correlate  flight  test  parameter  estimates 
with  wind  tunnel  data  and  analytic  results . 

2 .  The  desire  to  more  accurately  refine  parameter  estimates 
for  purposes  of  control  system  analysis  and  design. 

3 .  The  desire  to  achieve  an  accurate  aircraft  math  model  for 
use  in  high  fidelity  flight  simulators. 


An  early  and  continued  use  of  parameter  estimation,  as  stated 
above,  is  in  the  validation  of  wind  tunnel  and  analytic 
results.  However,  due  to  the  continuing  advances  in  aircraft 
design  and  performance  capabilities,  the  ability  to  accurately 
extrapolate  wind-tunnel  test  results  is  diminishing  and  a 
greater  emphasis  is  being  placed  on  flight  test  results. 
[Ref .l:p.2] 

Comprehensive  wind-tunnel  testing,  combined  with 
analytic  analysis,  can  give  reasonable  estimates  of  an 
aircraft's  aerodynamic  derivatives,  but  there  are  potential 
sources  for  inaccurate  predictions:  the  matching  of  "scaled" 
wind-tunnel   tests   with   expected   flight   conditions   is 


difficult,  with  Reynolds  number  differences  often  being  the 
standard  explanation  for  discrepancies .  Reliable  and  accurate 
dynamic  wind  tunnel  test  results  are  extremely  difficult  and 
expensive  to  accomplish.  Support  systems  (stings,  etc.)  have 
become  an  issue  as  to  the  extent  the  data,  especially  drag, 
are  affected.  [Ref.  2:p.3]  Additionally,  the  ever  present 
time  and  money  constraints  often  necessitate  shortcuts  in  not 
only  wind-tunnel  testing  but  in  flight  testing  as  well.  It 
seems  wise,  therefore,  to  use  flight-test  data,  at  the  very 
least,  as  a  verification  tool  of  aircraft  stability  and 
control  derivatives  for  even  the  most  simple  configurations. 

Currently  at  the  Naval  Postgraduate  School  (NPS) ,  in  the 
Department  of  Aeronautics  and  Astronautics,  research  is  being 
done  using  remotely  piloted  aircraft  as  research  testbeds . 
One  testbed  in  use  has  been  a  half-scale  Pioneer  Unmanned  Air 
Vehicle  (UAV) . 

The  full-scale  Pioneer  is  operational  in  the  U.S.  Navy  and 
Marine  Corps  and  saw  extensive  action  in  the  recent  war  with 
Iraq.  The  small  size  of  the  Pioneer  or  essentially  any 
tactical  UAV  allows  it  to  operate  close  to  and  in  some  cases 
behind  enemy  lines,  extending  the  "eyes"  of  battlefield 
commanders .  Its  missions  include  gun  fire  spotting,  real  time 
enemy  surveillance,  bomb  damage  assessment,  target  designation 
and  an  array  of  intelligence  collection  techniques. 

The  relatively  low  cost,  small  size,  reduced  risk  and 
inherent  flexibility  of  UAV's,  such  as  the  half-scale  Pioneer, 


have  allowed  the  department  to  become  actively  involved  in 
assessing  their  flight  characteristics.  The  Pacific  Missile 
Test  Center  (PMTC)  at  Pt .  Mugu  California  is  a  development  and 
testing  facility  for  the  U.S.  Navy.  Current  activity  at  PMTC 
includes  developmental  work  in  conjunction  with  the  Pioneer 
UAV.  In  development  by  the  Target  Simulation  Lab  at  PMTC  is 
a  flight  training  simulator  to  be  used  by  Pioneer  operators 
for  initial  training  and  proficiency  flights .  Results  from 
thesis  work  of  two  former  NPS  students,  USMC  Capts.  Daniel 
Lyons  and  Robert  Bray,  have  been  supplied  to  the  Lab  [Ref s .  3 
and  4]  .  These  results  comprised  various  aerodynamic 
parameters  obtained  from  two  different  approaches,  a  numerical 
method  (low-order  panel  technique)  and  wind  tunnel  tests  of  a 
0.4-scale  model.  The  aerodynamic  data  supplied  to  PMTC  are 
being  used  in  their  math  model  for  the  simulator. 

The  goal  of  the  present  study  is  to  incorporate  a  personal 
computer  (PC)  parameter  estimation  capability  into  the  ongoing 
NPS  flight  research.  This  application  will  give  the  flight 
test  program  an  added  dimension:  the  ability  to  compare  data 
from  wind  tunnel  and  analytic  analyses  with  flight  test 
results.  Additionally,  the  adapted  program  can  be 
incorporated  into  flight  test  and  dynamic  stability  and 
control  courses  as  a  valuable  teaching  aid. 

In  the  near  term,  interest  in  the  Pioneer  parameter 
estimation  results  has  been  expressed  by  PMTC  in  hopes  of 
achieving   a   more   realistic   training   simulator   for   the 


operators.  It  is  hoped  that  full  scale  time  history  data  can 
be  obtained  from  PMTC .  Future  work  in  this  area  includes 
completion  of  the  Pioneer  flight  research  and  comparison  of 
the  derivative  results  obtained  from  time  history  parameter 
estimation  with  those  obtained  in  References  3  and  4. 
Additionally,  other  UAV's  in  procurement  by  the  Department  of 
the  Defense  (DOD)  could  be  studied  at  NPS . 


II.    BACKGROUND 

A.   HISTORICAL  DEVELOPMENT 

The  history  of  flight  testing  would  in  itself  make  an 
alluring  and  fascinating  book;  this  cursory  summary  of 
parameter  estimation  and  flight  testing  reviews  but  a  small 
fraction  of  the  significant  events  in  the  history  of  flight 
testing.  The  majority  of  the  historical  content  was  found  in 
the  opening  remarks  given  by  Herman  A.  Rediess  [Ref .  5]  at  a 
1973  symposium  on  parameter  estimation  techniques. 

One  of  the  first  test  programs  to  obtain  quantitative 
measurements  of  aircraft  aerodynamic  characteristics  in  flight 
was  reported  by  Warner  and  Norton  in  1919  [Ref.  6]  .  Tests 
were  conducted  on  two  Curtiss  "Jenny"  JN-4H  biplanes  at 
Langley  Field,  Virginia.  Lift  and  drag  coefficients  were 
estimated  by  measuring  airspeed,  pitch  attitude  and  engine 
speed  in  flight  and  assuming  certain  engine  thrust 
characteristics.  This  specific  flight  test  was  a  meager 
beginning  for  in-flight  testing,  but  today  it  is  very  apparent 
that  in-flight  testing  is  a  vital  requirement.  A  1933  report 
by  Soule  and  Wheatly  [Ref.  7]  is  thought  to  be  the  first 
report  to  have  determined  and  compared  major  longitudinal 
stability  and  control  derivatives  from  flight  test  data  with 
results   acquired  through  theoretical   predictions.     The 


airplane  was  the  single  engine  Doyle  0-2.  The  analysis  used 
simplified  models,  solving  for  one  parameter  at  a  time  while 
assuming  values  for  the  other  parameters  based  upon  wind 
tunnel  data  or  other  flight  tests.  Early  in  the  1950' s  a 
major  advancement  in  parameter  estimation  was  achieved  by 
Shinbrot  [Ref.  8]  using  least  squares  curve  fits  between  the 
equations  of  motion  and  flight  data.  A  considerable  drawback 
at  that  time  was  the  extensive  calculations  required  for  this 
approach.  These  calculations,  of  course,  were  completed 
entirely  by  hand  as  the  digital  computer  was  not  yet 
available  as  an  engineering  tool. 

Significant  improvements  were  further  realized  in  the 
later  1950' s,  and  throughout  the  60' s  and  70' s,  due  to: 

1.  The  availability  of  the  digital  computer. 

2  .  The  progress  in  the  technical  disciplines  of  system 
identification  and  numerical  analysis. 

3.  The  availability  of  high  speed  automatic  data  acquisition 
systems . 

Again,  an  excellent  overview  of  the  evolution  of  parameter 
estimation  techniques  up  to  approximately  1970  is  contained  in 
Reference  5. 

There  are  numerous  methods  for  extracting  the  stability 
and  control  derivatives  from  flight-test  data  that  have  been 
developed  and  tested  since  the  early  50' s.  Each  starts  with 
equations  of  motion  and  essentially  attempts  to  curve  fit 
calculated  results  to  the  flight-test  data  by  adjusting  each 


of  the  derivatives  or  coefficients  in  the  math  model.  Some, 
but  certainly  not  all,  techniques  that  have  been  used  with 
success  include:  Ordinary  Least  Squares,  Weighted  Least 
Squares,  Deterministic  Least  Squares,  Maximum  Likelihood, 
Statistical  Linearized  Filter,  Extended  Kalman  Filter, 
frequency  domain  methods,  and  an  older  technique  used  in  the 
1950' s,  called  analog  matching.  This  last  technique  was  a 
manual  curve  fitting  method  using  an  analog  computer  in  which 
the  results  were  very  much  dependent  upon  the  skill  of  the 
operator . 

B.   MODERN  DEVELOPMENT 

Major  contributions  to  aircraft  parameter  estimation  since 
the  mid  1960's  have  been  made  at  two  NASA  facilities,  the  NASA 
Ames  Research  Center' s  Dryden  Flight  Research  Facility  and  the 
NASA  Langley  Research  Center.  The  parameter  estimation 
contributions  from  these  facilities  have  been  made  primarily 
through  the  work  of  Lawrence  Taylor,  Kenneth  Iliff,  Richard 
Maine  and  James  Murray. 

Taylor  and  Iliff  developed  a  Newton-Raphson  parameter 
estimation  program  in  1966  based  upon  the  theoretical  work  of 
Balakrishnan  [Ref.  9].  The  program  used  a  modified  Newton- 
Raphson  algorithm  to  effect  the  maximum  likelihood  technique 
for  estimating  stability  and  control  derivatives.  This 
program  underwent  a  gradual  evolution  during  its  application 
[Ref.  10].    The  outcome  in  1973  was  a  program  named  MMLE 


(Modified  Maximum  Likelihood  Estimator)  which  used  the  same 
basic  algorithm  (Newton-Raphson)  but  incorporated  features 
useful  for  processing  large  amounts  of  flight  data.  This 
program  was  widely  circulated  among  industry  and  government 
agencies  and  as  of  1980  had  been  used  to  analyze  over  35 
different  aircraft  at  the  NASA  Dryden  Flight  Research  Center 
alone.  [Ref.  ll:p.2] 

Development  of  MMLE3  (Modified  Maximum  Likelihood 
Estimation  program  version  3)  was  completed  by  Maine  and  Iliff 
in  1980  in  response  to  a  requirement  for  a  more  versatile 
parameter  estimation  program.  MMLE3  had  two  advances  over 
MMLE :  more  flexibility  in  defining  the  equations  of  motion 
(although  still  linear) ;  and  the  capability  for  estimation  in 
the  presence  of  state  noise,  also  called  process  or  input 
noise,  a  good  example  being  atmospheric  turbulence.  [Ref. 
11  :p.  2]   Further  details  on  MMLE3  will  be  addressed  later. 

In  1987,  a  new  parameter  estimation  program  to  accommodate 
nonlinear  models  was  developed  at  NASA  Dryden  by  Maine  and 
Murray  [Ref.  12].  This  parameter  estimation  program,  named 
pEst,  supports  nonlinear  models,  and  thus  aircraft  dynamic 
behavior  can  be  tested  at  extreme  flight  conditions  (high  AOA) 
and  for  unique  configurations  (oblique  wing,  etc.) . [Ref.  2] 

The  basic  concepts  of  aircraft  parameter  estimation 
techniques  have  remained  unchanged  for  over  two  decades  and 
are  shown  in  Figure  1  [Ref.  5: p. 14] .  These  concepts  include: 
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Figure  1.     Basic   Concepts   of   Contemporary   Parameter 
Estimation  Techniques 


(1)  the  mathematical  model,  (2)  the  data  acquisition  system, 
(3)   the  estimation  algorithm,   and   (4)   the  required  test 
inputs.   Each  of  these  elements  will  be  discussed  later  in 
more  detail . 

A  Maximum  Likelihood  (ML)  estimator  was  chosen  as  our 
parameter  estimator  for  the  following  reasons.  First,  this 
method  has  become  and  still  is  "accepted  as  the  standard 
approach  to  determining  aircraft  stability  and  control 
derivatives  from  flight  data"  [Ref . 13 :p . 558] .  Furthermore, 
the  flight  regimes  of  the  test  vehicles  at  NPS,  at  least 
initially,  are  expected  to  be  well  within  the  region  where  a 
linear  math  model  will  provide  accurate  parameter  estimations. 
Furthermore,  a  PC  compatible  ML  estimator  program  was  chosen 
because  of  the  PC's  flexibility  and  availability  at  NPS.  This 
combination  appeared  ideal  for  use  in  analyzing  the  ongoing 


NPS  flight  research  and  also  for  use  in  the  classroom  as  an 
enhancement  to  existing  teaching  aids. 
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III.   MATHEMATICAL  MODEL 

In  this  section,  the  linear  equations  of  motion  used  to 
describe  a  typical  manned  or  unmanned  aircraft  will  be 
developed.  Similar  developments  are  shown  in  most  aircraft 
dynamics  textbooks.  References  used  in  this  development  were 
Airplane  Flight  Dynamics  by  Roskam  [Ref.  14],  classroom  notes 
by  NPS  Professor  Louis  Schmidt  [Ref.  15],  the  textbook  Flight 
Stability  and  Automatic  Control  by  Nelson  [Ref.  16]  and 
Reference  2.  In  creating  simulated  data  a  gust  or  turbulence 
model  was  used  and  that  too  is  developed  in  this  section. 
Lastly,  a  discussion  of  the  observation  equation  corrections 
is  presented. 

A.   MODEL  EQUATION  DEVELOPMENT 

This  development  begins  with  Newton' s  linear  and  angular 
momentum  equations: 

F=JjL(mV)  (la) 

dt 

M=-ZL{H)  (lb) 

dt 

Where  the  force,  F,  is  the  sum  of  the  externally  applied 
forces  and  the  moment,  M,  is  the  sum  of  the  applied  moments 
about  the  center  of  gravity  (eg) .  The  use  of  non-rotating, 
earth  reference  coordinates   for  equations   la  and  lb  is 
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unwieldy  for  two  reasons.  First,  required  measurements  are 
predominately  made  in  the  rotating  body  axis  system;  and 
secondly,  but  of  more  significance,  the  inertia  matrix  or 
tensor  is  a  function  of  time  in  the  non-rotating  system. 
Therefore,  the  axis  system  chosen  is  the  standard  body  axis 
system,  shown  in  Figure  2  [Ref .  15] .  The  body  axis  system  is 
used  by  Iliff  and  Maine  in  Reference  2,  but  the  reader  should 
be  advised  that  the  equations  of  motion  are  also  at  times 
derived  using  the  stability  axis  system.  Reference  14  has  a 
good  description  of  the  differences  between  the  two  axis 
systems.  The  origin  is  positioned  at  the  vehicle's  eg  with 
the  X-direction  pointing  out  the  nose  of  the  aircraft,  the  Y- 
direction  out  the  starboard  wing  and  the  Z-direction  out  the 
bottom  of  the  aircraft .  Transforming  equations  la  and  lb  into 
the  rotating  body  axis  system  is  done  below: 

F=  3  (mv)  +wx  (mV)  (2a) 

ot 

M=-^-  (H)  +(oxH  (2b) 

ot 

where  the  angular  momentum,  H,  is  the  inner  dot  product  of  the 
aircraft  mass  moment  of  inertia  matrix,  [Imo]  ,  with  the  angular 
rotation  vector  co.  The  inertia  matrix  in  the  body  axis 
coordinate  system  is  not  a  function  of  time  as  it  is  in  the 
earth  reference.  The  above  equations  are  vector  equations  and 
can  be  written  into  scalar  components. 
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The  components  in  the  body  axis  for  the  CO  vector  are  p, 
q,  and  r  for  the  roll,  pitch  and  yaw  rates,  respectively.  The 
components  in  the  body  axis  for  V  are  u,  v,  and  w  for  the  X, 
Y,  and  Z  components  of  velocity,  also  shown  in  Figure  2. 

The  applied  forces  on  the  aircraft  can  be  broken  down  into 
aerodynamic,  gravitational  and  thrust  components.  (The  thrust 
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Figure  2 .   Body  Axis  System  and  Notation 

is  assumed  to  act  along  the  X-body  axis.) 

The  moments  can  be  broken  down  into  just  aerodynamic 
components,  since  the  thrust  is  assumed,  and  gravity  forces, 
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by  definition,  act  through  the  eg,  and  their  moment 
contributions  are  zero.  The  moment  components  are  then  due  to 
the  aerodynamic  forces  and  are  shown  in  Figure  2  as  L,  M,  and 
N. 

Since  the  gravity  force  components  in  the  body  axis  system 
depend  on  the  orientation  of  the  airplane  relative  to  earth- 
fixed  coordinates  (assuming  a  flat  non-rotating  earth) ,  it  is 
necessary  to  describe  the  orientation  of  the  aircraft  relative 
to  the  earth.  Euler  angles  are  introduced  to  accomplish  this 
transformation  between  coordinate  systems.  The  Euler  angles, 
*F,  0,  and  <£,  are  three  consecutive  rotation  angles  needed  for 
the  transformations  from  one  axis  system  to  the  other.  The 
angle  *¥  is  referred  to  as  the  yaw  angle.  The  angle  0  is  the 
pitch  angle.   The  angle  O  is  the  bank  or  roll  angle. 

The  aircraft  is  further  assumed  to  be  rigid;  that  is,  the 
mass  particles  remain  at  constant  distances  from  each  other. 
The  X-Z  plane  is  assumed  to  be  a  plane  of  symmetry  and  thus, 
the  products  of  inertia  Ixy  and  Iyz  are  zero.  In  this  model, 
the  rotating  engine  parts  and  sloshing  fuel  are  being  ignored, 
and  over  short  periods  of  time  when  data  are  collected,  the 
mass  of  the  aircraft  is  considered  to  be  constant. 

The  force  component  equations  that  are  derived  with  the 
above  assumptions  are: 
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Fx=m{u+qw+rv)  (3a) 


Fy=m{v+ru+pw)  (3b) 


Fz=m{w+pv+qu)  (3c) 

and  the  moment  component  equations  are: 

L=Mx=pIx-fIxz+qr(Iz-Iy)  -pqlxz  (4a) 

M=My=qIy+rp(Ix-Iz)  +  (p2-r2)  Ixz  (4b) 

N=Mz=-pIxz+fIz+pq(Iy-Ix)  +qrlxz  (4c) 

Where     the     left     hand     sides     of     the     above     force     and    moment 
component    equations    are: 

Fx=qsCx-mgsinQ  +  Thrust  (5a) 


Fy=qrsCy+/r?grsin<t>cos0  (5b) 


Fz=qsC z+mgcos$cosQ  (5c) 


Mx=gsjbCj  (6a) 
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My=qscCm  (6b) 


M2=qsbCn  (6c) 

where  q  is  the  dynamic  pressure  and  q  is  the  pitch  rate. 

The  vehicle  is  assumed  to  operate  at  small  side  slip 
angles  and  small  perturbations  around  a  steady  state 
condition.  The  perturbations  are  assumed  to  be  small  in  order 
that  the  sines  and  cosines  of  the  disturbance  angles  are 
approximately  the  angles  themselves  and  one  (1)  respectively, 
and  the  products  and  squares  of  the  products  are  negligible 
when  compared  to  the  quantities  themselves.  This 
approximation  is  termed  small  perturbation  theory  and  permits 
the  equations  of  motion  to  be  decoupled  and  linearized  into 
two  smaller  subsets:   Lateral-Directional  and  Longitudinal. 

It  is  often  times  more  convenient  to  have  the  equations  in 
terms  of  a  (angle  of  attack)  ,  |3  (sideslip  angle)  ,  and  V  than 
in  terms  of  u,  v,  and  w.  These  angles  are  usually  measured 
directly  vice  the  velocity  components.  In  the  transformation 
of  the  equations  of  motion  into  equations  with  a  and  ($,  the 
following  can  be  noted: 

a  =  tarrx(-)  (7a) 
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P  =  tan"1(-)  (7b) 

u 


and  if  the  angles  are  small,  then 


a*Z  (8a) 

u 


d«^  (8b) 

u 


p  =  ^  (8c) 


u 


p  =  ^  (8d) 

u 


These  equations  will  be  used  in  the  next  sections  to  construct 
the  basic  aircraft  model . 

B.   SIMPLIFIED  LONGITUDINAL  EQUATIONS 

The  longitudinal  set  of  equations  pertains  to  rotation  or 
moments  about  the  Y  axis  (4b  and  6b)  with  translation  or 
forces  along  the  X  and  Z  axis  (3a,  3c,  5a  and  5c) .  With  the 
substitution  of  8a  and  8b  from  above,  and  also  with  the  use  of 
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small   perturbation   theory,    the   following   simplified 
longitudinal  equations  are  formulated: 

d  =  ^C+g+-2  (9) 

mu  u 


fr-spq.  do) 

y 


6=g  (ID 

where  the  lift  is  approximately  parallel  to  the  Z-axis,  so 

CL*-CZ  (12a) 

^=CLaa+CLe6e  +  CLo  (12b) 


C  =C   a+C    -2£+C     6  +C  (12c) 


Equations  12a  and  12b  are  then  substituted  into  equation  9 
while  equation  12c  is  substituted  into  equation  10.  The 
result  expressed  in  state-space  format  and  using  dimensional 
derivatives  follows: 
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with  the  output  equation  being 
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In  equation  13a  the  Mq  term  includes  the  M^  term.  The  measured 
variables  are  V,  a,  q,  0,  A„  (normal  acceleration  in  "g's") 
and  6.  during  the  maneuvers.  The  magnitude  of  the  velocity, 
V,  is  approximately  equal  to  u  at  small  a  and  has  been 
substituted  for  u  in  the  above  state  space  equation.  The 
dimensional  derivatives  formulation  is  shown  in  the  list  of 
symbols . 

C.   SIMPLIFIED  LATERAL-DIRECTIONAL  EQUATIONS 

The  lateral-directional  equation  set  pertains  to  the 
rotation  about  the  X  and  Z  axes  (4a,  4c,  6a  and  6c)  and 
translation  along  the  Y  axis  (3b  and  5b)  .  These  equations  are 
used  to  derive  the  lateral-directional  state  space 
representation.  Small  perturbation  theory,  and  the  use  of 
equations  8c  and  8d,  are  used  to  obtain  the  following  lateral- 
directional  equations: 
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and 


Cy=CyJ+Cy**r  (14a) 


■y    ~ypr     ~yt[ 


$=qsCy+9$-r  (14b) 


mu        u 


1       iPr       jp27       j*  2^       Ifi«    a       ier    r 


plx-rlxz^qsbc^qr{l-lz)  +pqlxz  (15b) 


C  =c    6+C    ^  +  c    — +C    6  +C     6  (16a) 


-PIxz+rIz=qsbCn+pq(Ix-Iy)  -qrlxz  (16b) 


— ^L  =<j)=p+rtan0sin<|>+gtan0cos<|>  (17) 

at 


Substituting  equations  14a,  15a,  and  16a  into  14b,  15b  and  16b 
respectively,  and  again  expressing  these  equations  along  with 
equation  17  in  state-space  format  using  dimensional 
derivatives   the    following   is    obtained: 
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The  measured  variable  are  the  5a  and  5r  control  deflections, 
(3,  p,  r,  O,  Ay  and  V. 

The  two  state  space  representations  (longitudinal  and 
lateral-directional)  form  the  mathematical  model  used  with  the 
parameter  estimation  program  to  estimate  the  stability  and 
control  derivatives. 

D.   TURBULENCE  MODEL 

In  preparation  for  using  actual  time  history  data, 
simulated  longitudinal  and  lateral-directional  data  were 
created.  In  creating  the  simulated  data  it  was  desired  to 
match  the  expected  phenomena  or  real  world  effects  that  would 
be  encountered.  Thus  turbulence  (state  noise)  and  measurement 
noise  were  selectively  added  to  the  simulated  data  by  the 
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user.  This  section  describes  the  model  used  to  generate  the 
turbulence  or  state  noise. 

The  development  begins  with  the  application  of  the  Dryden 
gust  model  [Refs.  15  and  18]  .  The  turbulent  effect  can  be 
added  to  both  the  longitudinal  and  lateral-directional 
equations.  The  development  for  the  longitudinal  case  will  be 
shown  below  and  can  be  extrapolated  in  a  straightforward 
manner  for  the  lateral-directional  case. 

In  the  Laplace  domain  the  vertical  gust  velocity  transfer 
function  is 

wg{s)=JZ   ,S4*  T)(s)  (19) 

(s+X)2 


where 


*>=——  (20a) 


A,  =  ~-  (20b) 


*=2^  (20c) 
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and  L^,  is  the  scale  of  the  turbulence,  a„  is  the  rms  value  of 
the  turbulence  and  T|  is  a  zero  mean  white  noise  input .  The 
values  used  for  L^  and  <TW  equated  to  a  turbulence  level 
between  light  and  moderate,  and  can  be  adjusted  if  desired. 

Equation  21  was  obtained  by  transforming  equation  19  from 
the  Laplace  domain  into  the  time  domain  and  dividing  by  u: 


-^=ft  +2Adc+A2a  =J^  (iji|+1)  (2D 


u  9  9  9         u 


where  ag  is  the  a  perturbation  attributable  to  the  gust . 
Equation  21  is  converted  into  the  state  space  format  by 
letting 


Zx=a  (22a) 


Z2=d  -^t!  (22b) 

9      u 


Z^Z^&T)  (22c) 


Z7  =  -2XZ2-k2Z1  +  ^T)  (b-2X)  (22d) 


it  follows  that 
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This  result  can  be  combined  with  the  longitudinal  state  space 
equation  13  (Lateral-Directional  equation  18)  to  yield 
equation  24a: 


a 

e 


a 

V 

1 

0 

a 
V 

M* 

*** 

0 

tf« 

0 

1 

0 

0 

0 

0 

0 

0 

0   0  0  -\'< 


0 

V 

Q 

0 

e 

+ 

0 

z1 

1 

2k 

<Z2) 

V 


0 
0 


v^ 


0 

9 

V 

0 

0 

V 

0 

0 

Tl 

{R 

0 

UJ 

V 

(b-2k) 

0 

(24a) 


and  the  output  equation  24b: 
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Similarly,  the  state  space  equations  for  the  lateral- 
directional  model  with  turbulent  noise  can  be  developed  and 
are  shown  below: 
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E.   MEASUREMENT  CORRECTIONS 

The  observed  or  measured  data  must  be  corrected  for 
measurement  errors  caused  by  the  positioning  of  the  sensors. 
The  aircraft  equations  of  motion  were  developed  for  the 
aircraft  eg.  Therefore,  measurements  taken  by  sensors  not 
physically  located  at  the  eg  require  corrections  to  reference 
the  values  to  the  eg  location.   This  subsection  will  document 
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the  correction  equations  used  in  this  analysis.  The  simulated 
data  need  not  be  corrected  since  these  data  were  manufactured 
at  the  eg  of  the  aircraft.  However,  in  anticipation  of  actual 
aircraft  flight  test  parameter  estimation,  sensor  position 
corrections  were  implemented  into  the  programs  used  to  analyze 
genuine  flight  data. 
1 .   Longitudinal 

The  longitudinal  a  and  A„  data  are  capable  of 
corrections  for  sensor  position  displacement  from  the  eg.  The 
a  was  corrected  for  the  X-coordinate  probe  position  forward 
(  +  )  or  aft  (-)  of  the  eg  (Xap)  as  shown  below: 

a  =a  u.  +^a  (26a) 

eg       probe    T/ 

A  correction  for  the  upwash  angle  a*  at  the  probe  was  not 
taken  into  account  as  it  is  assumed  to  be  small  or  previously 
accounted  for  when  the  sensor  was  calibrated;  a  more  complete 
discussion  on  corrections  is  contained  in  Reference  2.  The 
normal  acceleration,  A„,  was  corrected  for  both  XM  (fwd  +)  and 
Zan  (down  +)  displacement  from  the  aircraft  eg  as  shown  below: 

A     =An     -^S-^Hq*  (26b) 

"cy    "ace    g  g 
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2 .   Lateral-Directional 

The  sideslip  angle,  |3,  and  the  lateral  acceleration, 
Ay,  were  corrected  for  their  sensor  positions  and  the 
correction  equations  are  shown  below: 

Pc-Pt*.-^  <27a) 


eg    rprobe         v 


and 


A     =A      -fltyi  +  ^vp  (27b) 


'  eg         ■>  ace 


The  values  for  Xbp,  Xay  and  Zay  are  again  defined  positive 
forward  and  down  analogous  to  the  body  coordinates . 


27 


IV.    MAXIMUM  LIKELIHOOD  PARAMETER  ESTIMATION 

A .   THEORY 

The  concept  of  parameter  estimation,  as  discussed  earlier, 
can  be  defined  quite  simply  in  general  terms.  The  system,  in 
this  case  a  UAV  or  some  other  testbed  whose  parameters  are  to 
be  estimated,  is  assumed  to  be  described  by  a  set  of  linear 
dynamic  equations,  a  mathematical  model,  which  was  defined 
earlier.  To  determine  the  values  for  the  unknown  parameters 
the  system  is  excited  by  an  input .  The  input  and  the  system' s 
actual  response  are  measured.  The  values  of  the  system 
unknowns  are  then  calculated  based  on  the  requirement  that  the 
model  response  to  the  input  match  the  actual  system  response. 
Complications  to  this  simplified  explanation  arise  when  the 
following  are  considered: 


1 .  Measurement  noise  -  perfect  measurements  are  unattainable 
with  any  sensor. 

2.  State  noise  -  the  aircraft  or  system  is  being  excited  by 
unmeasured  sources  such  as  atmospheric  turbulence. 

3.  Modeling  errors  -  exactly  describing  the  physical  system 
with  simple,  especially  linear,  dynamic  equations  is  very 
unlikely . 


In  fact,  if  the  above  complications  were  not  present,  the 
exact  values  of  the  unknown  parameters  could  be  found  and  what 
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is  termed  parameter  "identification"  vice  "estimation"  would 
be  our  accomplishment. 

The  common  approach  for  handling  the  modeling  error  is  to 
ignore  it  and  let  the  error  be  treated  as  measurement  or  state 
noise  or  both.  Iliff  and  Maine  state  that  "this  procedure  is 
not  rigorously  justifiable,  but  combined  with  a  carefully 
chosen  model,  it  is  probably  the  best  approach  available." 
[Ref.  19 :p.  2]  The  Modified  Maximum  Likelihood  Estimation 
algorithm  version  3  (MMLE3)  program  was  structured  to  take 
into  account  the  presence  of  state  and  measurement  noise. 
The  information  in  the  following  section  is  a  compilation  of 
information  on  MMLE  parameter  estimation  from  References  2  and 
19  through  22. 

B.   MODIFIED  MAXIMUM  LIKELIHOOD  ESTIMATION 

In  this  section  the  Modified  Maximum  Likelihood  Estimation 
algorithm  version  3  (MMLE3)  is  presented.  The  first  step  for 
parameter  estimation  then,  as  mentioned  above,  is  to  model  the 
system  accurately.  That  model  was  developed  in  the  previous 
chapter.  The  aircraft  equations  of  motion  define  the  system 
model  and  can  be  expressed  in  the  following  state  space  form: 

x(t0)=x0  (28) 


x{t)=Ax{t)+Bu(t)+Fw{t)+bx  (29a) 


29 


y(  t)  -Cx{t)  +Du(t)  +by  (2  9b) 


z{k)=y{k)  +Gv(k)  (29c) 

where  x(t)  is  the  state  vector  (x0  being  the  initial  state), 
u(t)  is  the  control  input  vector,  and  y(t)  is  the  prediction 
or  model  output  vector.  Matrices  A,  B,  C,  and  D  contain  the 
unknown  system  parameters,  which  in  this  case  are  the 
stability  and  control  derivatives.  Matrices  F  and  G  represent 
the  covariance  matrices  of  the  state  and  measurement  noise 
respectively.  The  measured  response  vector,  z (k) ,  is  sampled 
at  N  discrete  time  points  (k=l, . . .,N) . 

The  state  or  process  noise,  w(t),  is  assumed  to  be  zero- 
mean,  white  Gaussian  and  stationary.  The  measurement  noise, 
v(t),  is  assumed  to  be  zero-mean  Gaussian  noise  with  identity 
covariance . 

The  complete  unknown  parameter  vector  to  be  estimated  is 
then  given  by : 

{=(HT;\T;b?;b?)  (30) 

where  H  represents  the  unknown  parameters  in  the  matrices  A, 
B,  C,  and  D;  A.  represents  the  unknown  elements  of  the  F 
matrix;  and  bx  and  by  represent  the  unknown  biases  of  the  state 
and  model  output  equations  respectively.  The  [-]T  indicates 
the  transpose  of  a  matrix. 
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The  maximum  likelihood  estimates  are  obtained  by 
minimizing  the  negative  logarithm  of  the  likelihood  function. 
The  likelihood  function,  J,  is  a  function  of  the  difference 
between  the  measured  and  computed  time  histories .  The 
likelihood  function  is: 

J"(t»*)--!j;  [z(k)-y(k)]TR-1[z(k)-y(k))+^ln\R\         (31) 

2  Jc=l  2 


where  R  is  the  innovation  covariance  matrix.  The  innovations 
or  residuals  are  [z(k)-y(k)].  In  order  to  obtain  the 
predicted  state  variables  it  is  necessary  to  use  a  state 
estimator.  The  Kalman  filter,  which  is  an  optimal  linear 
state  estimator,  is  used  for  this  purpose.  The  Kalman  filter 
consists  of  a  prediction  step  and  a  correction  step  [Ref . 
20 :p.  12]  for  equations  2  9a  and  2  9b  and  is  shown  below: 

x{k+l)  =4>x(k)  +\\fBu(k)  +$bx  (32) 

y(k)=Cx(k)  +Du(k)+by  (33) 

x(k)=x(k)+K[z{k)-y(k)]  (34) 


U(/c)=-[u(*+l)+u(/c)]  (35) 

2 
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where  ~  (tilde)  and  (circumflex)  denote  the  predicted  and 
corrected  state  variables  respectively.  K  represents  the 
Kalman  filter  gain  matrix.  The  state  transition  matrix  is  $ 
and  its  integral  is  XP. 

1 .  Cost  Function  Minimization 

The  maximum  likelihood  estimates  for  the  unknowns  are 
found  by  minimizing  the  negative  logarithm  of  the  likelihood 
function,  J(^,R) .  The  negative  logarithm  of  the  likelihood 
function  is  often  called  or  referred  to  as  the  cost  function. 
This  minimization  is  done  by  using  a  modified  Newton-Raphson 
technique  which  iterates  on  the  vector  of  unknowns,  ^,  with 
each  iteration  providing  a  new  estimate  of  the  unknown  vector. 
These  new  estimates  update  the  math  model  coefficients, 
providing  a  new  calculated  response  and  a  new  response  error. 
This  iteration  process  is  continued  until  the  convergence 
criterion  is  satisfied. 

"The  maximum  likelihood  estimation  method  has  the 
desirable  characteristics  of  yielding  asymptotically  unbiased, 
consistent  and  efficient  estimates  [Ref.  19:p.3]." 

2 .  A-Priori  Weighting 

The  MMLE3  algorithm  allows  for  the  use  of  a  weighting 
function  to  account  for  prior  'engineering'  knowledge  of  the 
aircraft  parameters.  This  prior  knowledge  can  be  obtained 
from  other  test  cases,  wind  tunnel  measurements,  or  indepth 
analytic  analysis.  A  table  of  the  relative  importance  and  the 


32 


prediction  accuracy  of  stability  derivatives  using  theoretical 
methods  is  contained  in  Reference  14,  page  236  and  can  be  used 
as  a  guide,  if  desired,  to  weighting  the  initial  parameter 
estimates . 

The  a-priori  information  can  assist  the  algorithm  in 
converging,  but  caution  should  be  used,  as  the  weightings  can 
prejudice   the   answers   toward  the   analyst's   own   values 
[Ref .22:p.ST-8] . 

3 .   Estimate  Uncertainty 

The  use  of  the  Cramer-Rao  bounds  with  the  maximum 
likelihood  estimator  can  also  provide  a  measure  of  the 
relative  accuracy  of  the  estimates.  Each  parameter  bound 
gives  an  approximation  to  the  standard  deviation  of  the 
estimates.  It  is  important  to  recognize  that  these  bounds 
are,  in  fact,  the  lower  limits  for  the  standard  deviation, 
meaning  that  the  standard  deviation  value  is  at  least  as  large 
as  the  Cramer-Rao  bounds  [Ref.  21  :p.  12]  .  More  details  on  the 
Cramer-Rao  bound  is  contained  in  References  2,  20,  21  and  22. 
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The  State  Space  Identification  Toolbox  (SSIT)  for  the  386- 
MATLAB  personal  computer  program  implements  the  MMLE3 
algorithm.  MATLAB  is  a  registered  trademark  for  matrix 
oriented  software  distributed  under  license  agreement  by  The 
Mathworks,  Inc.  Reference  21  is  a  report  of  the  results  of  a 
study  comparing  the  mainframe  based  MMLE3  program  and  the 
MATLAB  SSIT  implementation  of  MMLE3 .  The  analysis  indicated 
that  the  PC  version  results  were  "generally  well  within  the 
uncertainty  levels  of  the  mainframe  parameter  estimates  [Ref . 
21:p.83]".  The  use  of  MATLAB  Software  will  be  discussed  in 
the  following  chapter. 
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V.    APPLICATION 

A.   DATA  ACQUISITION 

The  data  acquisition  system  is  an  important  part  of 
dynamic  stability  and  control  testing.  The  more  information 
known  on  the  details  of  the  entire  data  acquisition  system, 
the  greater  the  probability  that  the  test  results  can  be  more 
precise.  With  few  details  known  about  the  data,  often  times 
only  gross  characteristics  of  the  aircraft  can  be  determined. 
The  details  essential  for  a  complete  analysis  of  the  data 
should  include  how  the  data  were  filtered,  digitized,  time 
tagged,  transmitted  and  recorded.  The  complete  analysis  of 
the  data  acquisition  system  should  start  at  the  beginning  with 
the  sensors  and  continue  through  to  the  final  recorded  data 
product . 

Sensor  calibration  errors,  temperature  effects,  added 
noise  from  aircraft  vibration,  recorders  and  transmitters  are 
but  a  small  portion  of  the  circumstances  that  should  be 
reviewed.  Common  recording  systems,  their  advantages  and 
disadvantages,  plus  other  issues  relevant  to  the  entire  data 
acquisition  process  are  discussed  in  greater  detail  in 
Chapters  VII  and  VIII  of  Reference  2. 
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B.   INPUT  SELECTION 

The  selection  of  the  control  inputs  for  use  in  parameter 
estimation  must  take  into  account  the  pilot's  acceptability 
and  safety  of  flight  concerns  as  well  as  the  model  validity 
considerations . 

References  1,  2,  22  and  23  detail  various  methods  for  the 
input  design  employed  in  aircraft  parameter  estimation. 

One  specific  requirement  is  that  the  controls  applicable 
to  the  specific  model  need  to  be  exercised,  such  that  the 
aircraft  modes  are  excited.  Control  inputs  which  are  near  the 
frequency  of  the  excited  mode  usually  provide  the  best 
results.  This  is  because  at  these  modal  input  frequencies, 
the  largest  aircraft  response  for  a  given  input  usually  occurs 
and  provide  the  estimator  significant  data  as  compared  to  the 
noise  (gust  and  measurement  noise)  .  Judgment  must  be  used 
when  selecting  the  control  inputs  to  insure  flight  safety  and 
to  avoid  responses  that  exceed  any  preset  magnitude 
restrictions.  In  the  assumed  linear  model,  for  example,  as 
the  response  magnitudes  exceeded  the  small  perturbation 
assumptions  (10-15  degrees) ,  the  final  parameter  estimation 
results  can  be  expected  to  worsen.  Likewise,  the  signal-to- 
noise  ratio  of  the  data  improves  proportionally  with  the 
magnitude  of  the  response;  thus  there  is  a  trade-off  in  the 
development  of  the  aircraft  control  inputs.  Reference  23 
discusses  a  method  to  optimize  the  control  inputs  while 
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accounting   for   specific   restrictions   or   trade-offs   as 
mentioned  above. 

The  control  inputs  used  to  generate  the  simulated  data  for 
this  study  were  elevator,  rudder  and  aileron  ramped  doublets. 
This  type  of  control  input  was  selected  for  two  reasons. 
First,  this  input  is  truly  representative  of  actual  pilot 
inputs  (impulses  and  step  inputs  are  not  physically 
realizable) ,  and  second,  it  can  be  easily  adjusted  in 
magnitude  and  frequency  as  necessary  for  different  aircraft. 
The  previously  mentioned  references  provide  additional 
information  on  the  specifics  of  designing  control  inputs. 

C.   MATLAB 

MATLAB  is  a  commercially  available  software  package  for 
scientific  and  engineering  applications .  The  program 
integrates  numerical  analysis,  matrix  computation  and  graphics 
into  a  relatively  simple  environment  without  the  need  for 
traditional  programming  knowledge.  MATLAB  has  specialized 
toolboxes  for  added  capabilities .  In  this  study  the  Control 
Systems  Toolbox  and  the  State-Space  Identification  Toolbox 
(SSIT)  in  addition  to  the  main  386-MATLAB  program  were  used. 
The  SSIT  implements  the  MMLE3  algorithm,  discussed  previously 
in  Chapter  IV. 

Again,  the  use  of  the  MATLAB  based  program  is  desirable  at 
NPS  because  it  operates  in  a  familiar  PC  environment,  the 
importation  of  data  is  relatively  straightforward  and  easily 
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accomplished,  knowledge  of  a  formal  programming  language  is 
not  required,  and  the  plotting  and  hard  copy  functions  are 
easy  to  use  and  manipulate.  Furthermore,  during  the  required 
basic  dynamics  and  linear  systems  courses,  students  at  NPS 
use  MATLAB  extensively  as  an  instructional  and  problem-solving 
aid. 

1 .   M-f lies 

MATLAB  is  capable  of  executing  sequences  of  commands 
stored  in  files,  called  M-files  or  macros,  from  a  single-line 
command.  The  M-files  have  a  file  type  of  .m  and  consist  of  a 
sequence  of  normal  MATLAB  statements  that  can  include  the 
execution  of  other  M-files.  Major  benefits  of  the  M-files 
are:  repetitive  or  long  sequences  of  commands  can  be 
automated;  new  functions  can  be  created  by  the  user  for  a 
specific  need;  and  the  .m  files  are  ASCII  type  format  and 
easily  edited. 

The  following  sections  will  describe  the  M-files 
created  during  this  study  for  use  in  implementing  the  PC 
MATLAB  parameter  estimation  program.  It  is  noteworthy  that 
the  SSIT  is  itself  an  M-file  (mmle.m).  More  detailed 
information  concerning  MATLAB  is  contained  in  the  MATLAB 
user's  guides  References  22  and  24. 

a.  Simulated  Data 

Simulated  time  history  data  were  created  using 
MATLAB  M-files.   This  simulation  was  done  in  preparation  for 
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using  actual  aircraft  flight  data  in  the  parameter  estimation 
program.  The  aircraft  models  (longitudinal  and  lateral- 
directional)  used  to  create  the  data  were  previously 
discussed,  as  was  the  turbulence  model .  The  M-f iles  described 
in  the  following  sections  are  contained  in  Appendix  A. 
(1)    Longitudinal 

The  simulated  longitudinal  data  were  created 
using  the  M-file  LONGDAT.M  (see  Appendix  A) .  This  M-file  is 
extensively  commented  which  allows  the  user  to  understand  the 
program  and  change  or  adjust  certain  parameters  as  necessary, 
such  as  turbulence  level,  measurement  noise  level,  the 
elevator  input  amplitude,  and  period,  or  to  design  a 
completely  new  input . 

Aircraft  derivatives  and  other  physical  data 
are  necessary  to  create  the  simulated  data.  These  data  are 
stored  into  MATLAB  data  files  for  a  small  number  of  specific 
aircraft.  These  data  files  are  .mat  type  files.  The  user 
can  select  one  of  these  aircraft  or  input  the  data  for  an 
aircraft  of  his  or  her  choosing.  If  a  new  aircraft  is 
selected,  the  data  required  by  the  program  are  interactively 
requested  using  input  commands .  The  data  are  then  stored  and 
available  in  a  .mat  data  file  for  later  use.  The  storing  of 
these  data  into  accessible  files  eliminates  the  need  to 
reenter  the  data  every  time  additional  simulated  data  are 
desired  for  the  same  type  aircraft.    Aircraft  data  from 
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Reference  16  were  initially  entered  for  the  following 
aircraft:  NAVION,  A4D,  F104A,  JETSTAR,  and  B747 .  The  Pioneer 
UAV  data  were  also  entered  and  were  obtained  from  Reference  4. 
The  general  input  requirements  for  the  LONGDAT.M  macro  are: 


1.  The  aircraft  physical  data,  if  not  using  an  aircraft  with 
previously  saved  data 

2.  The  selection  of  either  adding  or  not  adding  state  and 
measurement  noise  to  the  data 

3.  The  flight  specifics:  velocity,  pressure  altitude,  and 
outside  air  temperature 


The  M-file  uses  a  ramped  elevator  doublet  for 
the  input  to  the  aircraft  math  model,  which  then  produces  the 
output.  The  simulation  can  be  done  either  with  or  without 
noise.  When  noise  is  selected,  in  addition  to  the  state 
turbulence  noise  added,  a  uniform  measurement  noise  is  also 
added  to  the  outputs.  The  measurement  noise  added  to  the 
outputs  a  and  0  is  zero  mean,  ±*i  degree  maximum  value,  while 
the  measurement  noise  added  to  the  output  q  is  zero  mean,  ±2.5 
degrees  per  second  maximum  value  and  the  measurement  noise 
added  to  the  normal  acceleration,  A^,  is  zero  mean,  ±.01  g 
maximum  value. 

The  final  output  time  histories  for  the  user 
include  5.,  a,  q,  0  and  A„  plots  displayed  on  the  PC  monitor, 
with  data  files  saved  containing  the  simulated  time  history 
data.   The  data  files  that  are  created  by  LONGDAT.M  are  then 
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available  for  the  parameter  estimation  M-files,  which  will  be 
discussed  later. 

(2)    Lateral-Directional 

The  simulated  lateral-directional  data  were 
created  using  the  macro  LATDIR.M.  This  file  is  similar  to 
LONGDAT.M  but  uses  the  lateral-directional  math  model.  The 
model  has  two  inputs,  5r  and  6a.  The  necessary  constants  and 
stability  and  control  derivatives  again  are  stored  into  a  .mat 
file  as  was  done  for  the  longitudinal  case.  The  noise  is 
essentially  the  same  as  in  the  longitudinal  model  except  that 
the  turbulent  gusts  are  caused  by  a  perturbation  in  the  side 
slip  angle,  |3g.  The  uniform  measurement  noise  that  is  added 
to  the  output  angles  ({3  and  <I>)  is  zero  mean,  ±H  degree  maximum 
value,  while  the  noise  added  to  the  output  angular  rates  (p 
and  r)  is  zero  mean,  ±2.5  degrees  per  second  maximum  value, 
and  the  noise  added  to  the  lateral  acceleration  is  the  same  as 
in  the  longitudinal  case. 

The  outputs  are  time  history  plots  and  data 
files  of  rudder  and  aileron  inputs,  5r  and  §«;  side  slip,  P; 
roll  rate,  p;  yaw  rate,  r;  roll  angle,  <X>;  and  lateral 
acceleration,  Ay.  The  data  files  that  are  created  by  LATDIR.M 
are  now  available  for  the  parameter  estimation  M-files,  which 
will  be  discussed  next. 
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b.      Parameter  Estimation 

Once  the  required  time  history  data  are  available, 
either  created  by  simulation  or  acquired  from  actual  flight 
tests,  the  parameter  estimation  algorithm  is  used  to  calculate 
the  'best'  estimate  of  the  stability  and  control  derivatives. 
The  process  of  executing  the  MATLAB  SSIT  parameter  estimation 
program  was  accomplished  with  four  basic  M-f iles .  The  four  M- 
files  were  developed  so  the  execution  of  the  SSIT  M-file 
(mmle.m)  would  appear  transparent  to  the  user.  Modifications 
to  the  basic  four  M-files  were  necessary  for  each  of  the 
following  cases: 

1 .  Simulated  longitudinal  data 

2 .  Simulated  lateral-directional  data 

3.  Actual  longitudinal  flight  data 

4 .  Actual  lateral-directional  flight  data 

These  M-files  are  included  in  Appendix  A.  The  four  basic  M- 
files  used  with  the  simulated  data  will  be  discussed  followed 
by  the  changes  needed  for  the  actual  flight  data  cases . 

The  four  basic  M-files  used  in  each  of  the  above 
four  cases  were  designed  to  simplify  the  execution  of  the  SSIT 
M-file  (mmle.m)  .  This  M-file  arrangement  is  shown  by  the 
block  diagram  in  Figure  3  below: 
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Figure  3  Block  Diagram  of  M-file  Arrangement 

The  four  cases  mentioned  above  each  have  the  four 
basic  M-files  shown  in  the  above  block  diagram.  The  above 
case  number  is  included  in  the  M-file  name  used  for  that  case 
to  distinguish  the  files  between  the  different  files.  Four 
files  for  each  of  the  four  cases  equate  to  16  different  M- 
files.  Each  set  of  files  is  a  variation  of  the  basic  four  M- 
f iles . 

The  only  M-file  the  user  initiates  is  the  macro 
named  NPSMMLE_  (the  is  where  the  case  number  from  above  is 
included  with  the  M-file  name)  shown  in  Figure  3.  The  other 
macros  are  "called"  in  a  manner  similar  to  that  for  a 
subroutine  call  in  FORTRAN  programming. 

NPSMMLE_  is  initiated  by  the  user  and  it  acts  as 
the  link  between  the  SSIT,  the  other  macros  and  the  user.  The 
one  exception  is  that  the  time  between  data  points,  dt,  must 
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be  edited  in  NPSP2SS_.M  for  cases  3  and  4.   The  functions  of 
NPSMMLE   include  the  following: 


1 .  It  initiates  the  MATLAB  diary  function  which  saves  a  copy 
of  the  output  (no  plots)  for  subsequent  printing 

2.  It  "calls"  the  NPSINIT_  macro  (discussed  later) 

3.  It  arranges  the  time  history  data  into  the  required 
format 

4 .  It  establishes  and  makes  known  the  function  file  to  the 
SSIT,  which  converts  the  parameter  vector  into  the  user 
defined  state  space  description  of  the  model 

5.  It  defines  additional  information  to  the  SSIT  which  is 
discussed  in  detail  in  References  21  and  22 

6.  It  "calls"  the  SSIT  parameter  estimation  program  (mmle.m) 

7 .  It  formats  and  displays  the  final  numerical  results  to 
the  monitor 

8.  It  "calls"  the  MLEPLOT_  macro  which  graphically  displays 
the  results 


The  NPSINIT_  macro  is  "called"  by  the  NPSMMLE 
macro  and  performs  the  following  functions: 


1.  It  loads  the  time  history  data  file  supplied  by  the  user. 
This  data  file  is  either  obtained  by  simulation  or  from 
actual  flights. 

2.  It  requests  the  time  interval  between  data  points. 

3.  It  requests  the  vehicle's  physical  attributes,  initial 
parameter  estimates  and  the  flight  conditions. 


All  of  the  above  information  is  saved  and  then  made  available 
for  the  SSIT  when  it  is  reloaded  into  the  MATLAB  working 
environment  by  NPSMMLE  . 
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The  NPSP2SS_  macro  is  a  function  file  used  by  SSIT 
during  the  parameter  estimation  process.  The  file  converts 
the  parameter  vector  into  the  user  defined  stated  space 
description  of  the  system.  Details  of  the  P2SS  M-file  are 
contained  in  References  21  and  22. 

After  the  SSIT  has  been  "called"  and  successfully 
run,  the  numerical  results  are  output  to  the  analyst  by  way  of 
the  PC  monitor.  The  MLEPLOT_  macro  is  then  initiated  and  its 
function  is  to  display  the  results  graphically  to  the  user. 
These  plots  are  shown  on  the  monitor  and  saved  as  MATLAB  .MET 
files.  These  .MET  files  are  high  resolution  graphics  files 
that  may  be  used  later  for  printing  graphics  hard  copies . 

The  macro  file  versions  numbered  1  or  3  are  for 
use  with  the  longitudinal  data,  simulated  and  actual  data 
respectively.  The  versions  numbered  2  or  4  are  for  use  with 
the  lateral-directional  data,  simulated  and  actual 
respectively . 

The  differences  between  the  versions  1  and  3  and 
2  and  4  are  relatively  small.  The  biggest  differences  are 
that  when  the  actual  data  (3  and  4)  are  used,  the  user  is 
asked  for  sensor  position  corrections  and  queried  whether  an 
a-priori  weighing  vector  is  to  be  used  during  the  parameter 
estimation  process  by  the  SSIT.  If  selected,  the  weighting 
input  is  a  vector  of  the  variances  for  each  of  the  parameters 
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to  be  estimated.  Therefore,  the  higher  the  number,  the  less 
the  weighting  afforded  that  parameter  and  vice  versa.  These 
M-files  are  included  in  Appendix  A. 
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VT.    RESULTS  AND  DISCUSSION 

A.   SIMULATED  DATA 
1 .   Application 

The  MATLAB  M-files  LONGDAT.M  and  LATDIR.M  were  used  to 
develop  the  required  simulated  data  representing  time  history 
responses.  These  data  included  both  state  and  measurement 
noise.  Qualitatively,  these  data,  plotted  in  Appendix  B  and 
discussed  below,  compare  favorably  with  actual  time  history 
data  shown  in  many  references.  Thus,  the  model  chosen  for  the 
simulation  was  assumed  as  an  adequate  mathematical 
representation  of  the  physical  system  within  the  region  of 
applicability,  this  region  being  the  area  or  flight  regime 
satisfying  the  restrictions  placed  upon  the  math  model  during 
development.  Additionally,  with  the  simulation,  there  is  the 
capability  for  the  user  to  modify  the  control  inputs  and  noise 
magnitudes  thereby  enhancing  the  use  of  these  programs  as  an 
instructional  aid.  Various  aircraft  types  with  differing 
inputs  and  atmospheric  conditions  potentially  could  be 
examined. 

The  parameter  estimation  of  the  simulated  data  was 
quite  accurate  with  some  exceptions.  These  exceptions  will  be 
discussed  later.  The  accurate  parameter  estimates,  for  the 
given  model,  validate  the  MMLE  methodology  in  the  presence  of 
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both  state  and  measurement  noise.    The  results  show  the 
relative  insensitivity  of  the  MMLE  method  (implemented  by  the 
MATLAB  SSIT)   to  noise  when  the  modeling  and  the  excitations 
(inputs)  are  chosen  with  care. 

Longitudinal  and  lateral-directional  parameter 
estimation  examples,  plots  and  quantitative  output  using  the 
simulated  data  are  shown  in  Appendix  B.  These  simulated 
examples  are  for  the  A-4D  and  Navion  aircraft  and  the  PIONEER 
UAV.  The  SSIT  quantitative  results  are  tabulated  with  the 
initial  input  parameter  estimates  and  the  "truth"  or 
underlying  parameters  (parameters  used  to  generate  the  data) 
for  comparison.  The  tabulated  SSIT  quantitative  results  are 
presented  in  column  format  with  the  column  headings  as 
follows:  pid,  parameter  id  number;  p (pid) ,  final  parameter 
estimate;  pref,  initial  input  reference  parameter;  cramer, 
Cramer-Rao  bounds;  2f cramer,  two  times  a  corrected,  filtered 
Cramer-Rao  bound;  insens,  insensitivity,  the  change  in  that 
parameter  required  to  move  from  the  maximum  likelihood  value 
to  the  edge  of  the  confidence  ellipsoid.  For  a  single 
parameter  model  the  insens  value  is  the  Cramer-Rao  bound. 
2 .   Longitudinal 

a. .       A-4D 

The  SSIT  parameter  results  for  the  A-4D,  due  to 
the  elevator  control  input  shown  in  Figure  B-l,  are  presented 
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in  Figures  B-2  through  B-5 .  These  are  the  input  and 
subsequent  response  of  the  A-4D .  The  a,  q,  0,  and 
acceleration  responses  correlate  well  with  the  truth  data  used 
to  generate  the  simulated  data.  All  estimated  parameter 
values  for  CL<X,  C^,  C^,  Cl5.,  and  C^.  shown  in  Appendix  B  are 
within  the  2fcramer  bound.  This  bound  represents  the  95 
percent  probability  ellipsoid  for  the  parameters.  Close 
inspection  of  the  estimates  show  that  the  greatest  deviations 
are  for  C^,  and  C^  and  both  are  within  four  percent  of  the 
actual  underlying  values  used  to  create  the  simulation. 

b.  Navion 

The  SSIT  results  for  the  Navion  are  shown  in  the 
response  plots,  Figure  B-7  through  B-10,  due  to  the  elevator 
control  doublet  input  shown  in  Figure  B-6.  Again,  as  was  seen 
in  the  previous  A-4D  results,  the  correlation  of  the  estimated 
aircraft  response  plots  with  the  underlying  truth  derivative 
responses  are  good.  The  largest  errors  are  in  the  C„q  and  Cl6<> 
derivatives,  and  they  are  12.4  and  7  percent  respectively.  An 
accurate  prediction  of  C^  using  theoretical  methods  is 
difficult  to  achieve;  Roskam  [Ref .  14]  notes  an  acceptable 
estimated  prediction  accuracy  for  this  derivative  of  20 
percent . 

c.  UAV 

The  results  for  the  A-4D  and  Navion  are  more 
accurate  than  the  UAV  results  which  are  shown  in  Figures  B-ll 
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through  B-15 .  The  UAV  results  demonstrate  the  importance  of 
designing  adequate  inputs  for  each  specific  vehicle.  The 
input  must  sufficiently  excite  the  vehicle's  dynamic  response. 
Identical  elevator  input  doublets  were  used  for  all  three 
vehicles;  the  UAV  elevator  input  is  shown  in  Figure  B-ll.  The 
input  period  is  2  seconds  with  a  maximum  amplitude  of 
approximately  4  degrees.  The  small  UAV  response  due  to  the 
small  elevator  control  power  (C^,)  might  be  overlooked  if  the 
scaling  on  the  response  plots  were  not  closely  observed.  The 
responses  are  less  than  a  half  to  a  third  those  of  the  A-4D 
and  Navion.  These  small  responses  equate  to  a  lower  signal- 
to-noise  ratio  for  the  parameter  estimator.  The  response  due 
to  the  elevator  input  was  not  much  more  significant  than  the 
response  perceived  by  the  estimator  from  the  state  noise  (the 
gust)  in  combination  with  the  presence  of  the  measurement 
noise.  Thus,  the  parameters  with  the  exception  of  CLa  (2.4 
percent)  were  all  in  error  greater  than  30  percent,  and  C„q  was 
76  percent  in  error  from  its  underlying  truth  value.  Caution 
must  be  exercised  in  choosing  a  proper  excitation  tailored  for 
the  particular  vehicle's  response. 

3 .   Lateral-Directional 

The  lateral-directional  parameter  estimates  provided 
by  the  SSIT  for  all  three  simulated  aircraft  are  also  shown  in 
Appendix  B . 
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A-4D 


The  A-4D  lateral-directional  inputs  are  shown  in 
Figure  B-16  and  consist  of  a  rudder  doublet  followed  by  an 
aileron  doublet.  The  aircraft  responses  of  6,  r,  <l>,  p,  and 
lateral  acceleration  due  to  the  aileron  and  rudder  inputs  are 
shown  in  Figures  B-17  through  B-21.  As  can  be  seen  in  these 
plots,  the  correlation  between  the  measured  and  estimated 
responses  is  good  and  the  parameter  estimates  which  are  also 
shown  in  Appendix  B  were  accurately  estimated.  Of  the  12 
parameters  all  but  three  were  inside  the  estimator  2fcramer 
bound.  The  three  derivatives  (i.e.  C^.,  Cl5a  and  CN8a)  ,  although 
not  inside  the  2fcramer  bound,  were  very  close  to  the  bound, 
and  are  within  18,  three  and  eight  percent  of  the  truth  values 
respectively.  Roskam  [Ref.  14  p. 236]  indicates  that  the 
theoretical  accuracy  for  values  in  estimating  CNr  is 
approximately  25  percent  using  analytical  methods .  It  is  felt 
that  improvements  in  the  estimation  of  these  derivatives  could 
be  achieved  by  investigating  different  inputs. 
b.       N avion 

The  Navion  results  are  also  shown  in  Appendix  B, 
Figures  B-23  through  B-27  are  the  response  plots  due  to  the 
rudder  and  aileron  inputs  as  described  by  Figure  B-22 .  The 
inputs  are  identical  to  the  rudder  and  aileron  inputs  used  for 
the  A-4D.    The  Navion  aircraft  responses  6,  r,  <t>,  p,  and 
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lateral  acceleration  to  the  inputs  (Figures  B-23  through  B-27) 
also  show  good  correlation  and  indicate  that  the  parameter 
estimates  were  estimated  accurately.  A  vertical  shift  in  the 
estimated  plot  from  the  true  response  sometimes  occurs  and  can 
be  misleading.  An  example  of  this  vertical  shift  is  shown  in 
Figure  B-26  by  the  Navion  roll  response.  This  misleading 
vertical  shift  is  caused  by  noise  on  the  first  data  point. 
The  estimated  plots  were  originated  by  MATLAB  at  the  first 
data  point  and  any  noise  in  that  point  causes  a  vertical  shift 
of  the  estimated  curve.  A  skilled  analyst  can  adjust  the 
first  point  to  the  known  initial  flight  condition  and 
eliminate  the  misleading  vertical  offset.  The  numeric  values 
for  the  Navion  parameter  estimates  are  also  shown  in  Appendix 
B  and  were  accurately  estimated.  Of  the  12  parameters  all 
were  inside  the  2fcramer  bound. 
c.   UAV 

The  parameter  estimation  results  for  the  UAV  due 
to  the  rudder  and  aileron  doublets  were  also  quite  accurate. 
The  aircraft  inputs  and  response  plots  are  shown  in  Figures  B- 
28  through  B-33 .  The  identical  rudder  and  aileron  inputs  used 
for  the  A-4D  and  Navion  were  used  for  the  UAV  and  are  shown  in 
Figure  B-28.  The  misleading  vertical  shift  is  more  prominent 
in  the  UAV  responses,  especially  in  Figures  B-30  and  B-33,  the 
yaw  rate  and  lateral  acceleration  plots.  Again,  the  responses 
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correlate  well   (except  for  the  vertical  shift)   and  the 
estimated  parameters  are  all  within  the  SSIT  2fcramer  bounds. 
4 .   Problems 

Significant  problems  were  encountered  when  trying  to 
employ  the  SSIT  on  a  very  lightly  damped  or  divergent  system. 
This  case  was  experienced  with  the  F104A  aircraft.  During 
these  conditions  the  SSIT  mmle.m  program  would  halt  prior  to 
completion  due  to  an  error.  The  error  displayed  to  the 
analyst  was  always  a  matrix  singularity  problem  encountered 
during  the  SSIT  computations. 

Another  rarely  encountered  problem  was  the  case  where 
the  calculated  estimates  were  significantly  in  error  from  the 
actual  parameters.  This  case  was  very  obvious  to  the  user, 
especially  when  the  estimated  plot  was  compared  with  the 
measured  data.  These  significant  parameter  errors  (1  to  2 
orders  of  magnitude)  caused  the  estimated  plot  to  rapidly 
diverge  from  the  measured  data.  It  is  believed  that  the 
estimation  program  was  converging  upon  a  local  minimum  instead 
of  the  global  minimum  of  the  cost  function. 

B.   ACTUAL  FLIGHT  DATA 
1 .   Application 

The  two  actual  flight  tests  analyzed  were  for  the 
longitudinal  cases  of  the  F-14A  and  the  T-37  aircraft.  The  F- 
14A  data  were  acquired  from  the  Naval  Air  Test  Center  and  the 
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T-37  data  from  NASA-Ames  Dryden  Flight  Research  Center.   The 
quantitative  results  and  plots  are  shown  in  Appendix  B. 

a.  F-14A 

The  F-14A  elevator  input  is  shown  in  Figure  B-34 
with  the  a,  q,  0  and  acceleration  responses  shown  in  Figures 
B-35  through  B-38 .  The  estimated  response  period  appears  to 
correlate  well  but  the  magnitude  of  the  estimated  responses 
are  less  than  the  actual  measured  responses  in  all  cases. 
This  result  could  be  due  to  the  accuracy  of  the  aircraft 
physical  characteristics  used  in  the  aircraft  model  or  to  data 
sensor  location  inaccuracies.  These  possible  problems  and  the 
predicted  derivatives  will  be  discussed  later.  The  estimation 
does,  however,  give  a  general  representation  of  the  aircraft 
longitudinal  characteristics. 

b.  T-37 

The  T-37  elevator  input  is  shown  in  Figure  B-39. 
The  aircraft  response  is  shown  in  Figures  B-40  through  B-43. 
The  response  plots  are  very  accurate  with  the  exception  of  a 
disparity  in  the  0  response  (Figure  B-42)  after  approximately 
the  5  second  point.  It  is  not  known  what  caused  this 
perturbation  since  the  q  (pitch  rate)  response  correlates  well 
and  is  the  time  derivative  of  0  response.  A  possible  cause 
could  be  attributed  to  a  data  acquisition  problem,  and  quite 
possibly  the  cause  will  never  be  known. 
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2 .   Problems 

The  problems  encountered  using  actual  flight  data 
involved  obtaining  accurate  aircraft  physical  characteristics 
and  arranging  the  data  into  the  desired  format.  Since  the 
flight  tests  were  not  specifically  designed  and  conducted  for 
this  study,  much  of  the  required  physical  data  were  not 
readily  available  and  were  estimated.  The  use  of  estimates 
was  especially  true  for  the  F-14A  data.  Additionally,  with  no 
truth  data  and  so  many  variables  on  the  F-14A  such  as  wing 
sweep  angle,  independent  left  and  right  horizontal  control 
surfaces  and  unknown  external  loadings,  the  accuracy  of  the 
model  used  and  the  results  are  at  best  an  approximation. 

In  both  cases  the  output  parameter  estimates  were 
obtained  without  the  use  of  the  program's  weighting 
capability.  The  parameter  estimates  appear  to  be  reasonable 
with  the  exception  of  the  CLa  values,  which  are  inclined 
towards  the  higher  side  of  discretion.  However,  with  only  one 
set  of  data  (one  flight  maneuver)  for  each  airplane  analyzed, 
no  significant  numerical  conclusions  can  be  justifiably 
deducted.  Continued  experience  with  actual  flight  data  will 
provide  a  database  upon  which  decisions  of  proper  weighting 
values  can  be  determined. 
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VII.     CONCLUSIONS  AND  RECOMMENDATIONS 

A.   Conclusions 

The  following  conclusions  were  deduced  from  this  study: 


1 .  The  simulated  time  history  data  generation  and  its  use  in 
the  parameter  estimation  program  worked  satisfactorily.  The 
accurate  parameter  estimates  validated  the  MMLE  method  as 
implemented  by  MATLAB  and  its  relative  insensitivity  to 
state  and  measurement  noise  if  the  model  and  inputs  were 
carefully  selected. 

2.  It  is  thought  that  significant  benefits  can  be  achieved 
with  the  simulated  portion  of  the  study  as  a  classroom 
instructional  aid  in  the  Flight  Test  and  Flight  Dynamics 
courses  at  NPS. 

3 .  The  actual  flight  test  data  appeared  to  produce 
acceptable  results.  However,  many  details  were  not  known 
concerning  the  data  and  aircraft  characteristics.  These 
unknown  details  concerned  such  items  as  data  filtering, 
sensor  position,  accurate  aircraft  physical  characteristics, 
and  time  lags.  Direct  involvement  with  the  flight  tests, 
although  not  essential,  would  have  significant  benefits. 
These  benefits  would  include  knowledge  of  the  above  missing 
details,  easier  data  acquisition  into  the  desired  data 
formats  and  hopefully  a  more  in-depth  and  accurate  analysis. 


B.   Recommendations 

The  following  recommendations  are  offered  based  upon  the 
above  conclusions: 


1 .  Incorporate  this  study  into  the  Flight  Test  and  Flight 
Dynamics  courses  at  NPS  as  instructional  demonstrations. 

2.  Continue  to  develop  an  on-site  data  reduction  capability 
to  enhance  the  flight  test  research  being  conducted  with 
UAV  s  at  NPS .  This  development  will  enable  the  school  to 
have  the  capability  of  comparing  computational  computer 
studies  and  wind-tunnel  studies  with  flight -test  results. 
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3.  Investigate  the  possibility  of  incorporating  the  pEst 
non-linear  parameter  estimation  program  (Reference  12)  into 
the  flight  test  research  program  to  cope  with  any  future 
requirements.  These  applications  could  include  helicopter 
dynamics,  high-a  flight  regimes  or  a  host  of  other  areas  not 
particularly  well  suited  to  linear  modeling. 

4 .  Investigate  developing  a  neural  network  parameter 
estimation  capability.  This  capability  could  be  used  for 
development  of  a  real-time  reconf igurable  flight  control 
system  to  improve  aircraft  survivability.  These  development 
areas  tie  parameter  estimation  into  two  additional  research 
disciplines  of  interest  at  NPS,  Neural  Networks  and  Aircraft 
Combat  Survivability. 

5.  Investigate  assimilating  the  GAT-1  training  device  and 
the  results  from  parameter  estimation  tests  into  a 
reconf igurable  simulator  for  instructional  purposes  in 
advanced  Flight  Dynamics,  Control  and  Avionics  courses. 

6.  Investigate  the  feasibility  of  using  the  MATLAB  SSIT  to 
accurately  determine  ship  dynamics  and  perhaps  find  better 
ways  of  controlling  unwanted  ship  motion.  The  improvement 
in  ship  dynamics  could  lead  to  improved  helicopter  landing 
conditions  during  rough  sea  states. 

7 .  Finally,  continue  the  NPS  flight  test  research  program  on 
Department  of  Defense  UAV  s  such  as  the  Pioneer  and  Exdrone, 
not  only  to  assist  in  their  evaluation  but  to  stimulate  the 
students'  interest  with  relevant  and  available  military 
research  topics . 
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APPENDIX  A  —  M-FILES 


LONGDAT.M 


clear; 

%   MACRO  FILE  NAME 

% 

r 


>=======  LONGDAT.M  =======< 

Date:  1  Feb  92 


disp 
disp 
disp 
disp 
disp 
disp 


disp 
disp 
disp 
disp 
disp 
disp 
disp 
disp 


ans) 
'  ') 


ans) 
'  ') 


MACRO  TO  GENERATE  SIMULATED  LONG.  DATA 
USING  ELEVATOR  DOUBLET  WITH  OR  W/O  NOISE 


') 


') 


GET  AIRCRAFT  TYPE  TO  BE  USED  

'  AIRCRAFT  TYPES  AVAILABLE  ' ) 

'  ') 

'NAVION    A4D    F104A     JETSTAR     B747     UAV   OTHER 

'  ') 

SELECT  OTHER  TO  INPUT  DATA  FOR  USER  DEFINED  AIRCRAFT 


') 


') 


') 


'NOTE  ========>  PROGRAM  IS  CASE  SENSITIVE  <====== 

'  ') 

typac=input ('TYPE  IN  DESIRED  A/C  FROM  THE  ABOVE  LIST.    ','s'); 

%  DETERMINE  IF  NOISE  IS  WANTED 

dispC  ') 

sysn=input (' INPUT  A  1  TO  INCLUDE  NOISE  AND  A  ZERO  FOR  NO  NOISE. 

ndp=201;%        NUMBER  OF  DATA  POINTS  -10  SEC 

dt=.05;%         TIME  STEP  FOR  THE  DATA 

amp=.07/%        AMPLITUDE  OF  DOUBLET  (RADIANS) 

period=l;%       PERIOD  OF  DOUBLET  IN  SEC=  PERIOD  +  1  SEC 

t=[0 :ndp-l] *dt;  %  TIME  VECTOR 

simdata=zeros (ndp,5) /         %  SETUP  DATA  MATRIX  ALL  ZEROS 

% GENERATE  THE  INPUT  DOUBLET 

dslope= (4*amp*dt) /period; 

dl=  (0  :  -l*dslope  :  -l*ainp)  ;dla=-amp*ones  (1:10)  ; 

d2= (-l*amp+dslope : dslope : amp) ; d2a=-l*dla; 

d3= (amp-dslope : -l*dslope : 0) ; 

simdata ( : , l)=[dl  dla  d2  d2a  d3  zeros (1, ndp- (period/dt) -21) ]' ; 

% GENERATE  THE  STABILITY  AND  CONTROL  MATRICIES 

dispC  ' ) 

vtrue=input (' INPUT  AIRCRAFT  TRUE  VELOCITY       ft/sec 

altft=input (' INPUT  AIRCRAFT  ALTITUDE  IN         ft 

oat =input (' INPUT  THE  OUTSIDE  AIR  TEMPERATURE    Deg  F 

if  (strcmpC  OTHER'  ,typac)>0)  ; 

typac=input (' INPUT  THE  A/C  TYPE  ....<  6  characters 
sref=input (' INPUT  THE  AIRCRAFT  REFERENCE  AREA    ft~2 
gw=input (' INPUT  THE  AIRCRAFT  GROSS  WEIGHT        lbs 


') 
') 
') 


s') 
') 
') 


)/ 
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iyy=input ('AIRCRAFT  Iyy  MOMENT  OF  INERTIA    slug-ftA2     '); 

cbar=input ('AIRCRAFT  MEAN  CHORD  LENGTH  ft       '); 

CLa=input (' INPUT  DERIVATIVE  CL_a      1/RAD 

CMa=input (' INPUT  DERIVATIVE  CM_a      1/RAD 

CMq=input (' INPUT  DERIVATIVE  CM_q      1/RAD 

CLde=input (' INPUT  DERIVATIVE  CL_de    1/RAD 

CMde=input (' INPUT  DERIVATIVE  CM_de    1/RAD 

ans=['save   ' , typac, ' .mat ; ' ] ; 

eval(ans);   %  SAVE  THE  NEW  A/C  DATA  IN  A  .MAT  FILE 

truth=[CLa  CMa  CMq  CLde  CMde]; 
else 

ans=['load  f , typac, ' .mat; f ] ; 

eval (ans) ; 

truth=[CLa  CMa  CMq  CLde  CMde]; 
end 

%    CALCULATE  DENSITY,  DYNAMIC  PRESSURE,  AND  CONSTANTS 
rho=.0023769*exp( (-1*32 . 17*altft) / (1716* (oat+460) ) ) ; 
qbar= . 5*rho* (vtrueA2) ; 

constl=-l* (qbar*sref ) / (gw*vtrue/32 . 17) ; 
const2=qbar*sref *cbar/iyy; 
const3=const2*cbar/ (2*vtrue) ; 
const 5=qbar*sref /gw; 
lw=1750;  %  SCALE  OF  TURBULENCE 

if  altft<1750 
lw=altft; 

end 
sigw=08;  %  RMS  VALUE  OF  TURBULENCE  IN  FT /SEC 

lamda= (vtrue/lw) ;k= (3*sigwA2) *vtrue/ (pi*lw) ;beta=vtrue/ (sqrt (3) *lw) 


% 


a= [constl*CLa 

1 

0 

constl*CLa 

0 

const2*CMa 

cons 

it  3*  CMq 

0 

const2*CMa 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

1 

0 
% 
b=[constl*CLde 

0 

0 

-(lamdaA2) 

-2*lamda] 

0 

0; 

const2*CMde 

0 

0, 

0 

0 

0, 

0 

sqrt (k) /vtrue 

0, 

0 
% 
c=[l 

(s 

sqrt (k) /vtrue) 

* (beta-2*lamda) 

0], 

0 

0 

0 

0; 

0 

1 

0 

0 

0; 

0 

0 

1 

0 

0; 

const5*CLa 
% 
d=[0 

0 

0 

0 

0]; 

0 

0; 

0 

0 

0; 

0 

0 

0; 
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const5*CLde    0        1]  ; 

% 

%SIMDATA(  :  ,2)=A0A  (RAD)    SIMDATA( :  ,  3)  =PITCH  RATE  (RAD/SEC) 

%SIMDATA( : , 4)=THETA  (RAD)  S IMDATA  (  :  ,  5 )  =NORMAL  ACC  (G) 

%     RAND  NUMBER  GENERATOR    MEAN=0      VARIANCE=1 

rand ( ' normal ' ) ; 
[phi, gam] =c2d (a,b, dt) ;u= [simdata ( : , 1) , sysn*rand (ndp, 1) , ones (ndp, 1) ] 

r 

[ynoise, xnoise] =dlsim (phi, gam, c, d, u) ; 

%  STATENOISE     +    MEASUREMENT  NOISE 

rand ( ' uniform' ) ; 

simdata ( : , 2 : 3 : 4 : 5) =ynoise (:, 1:2:3:4)  ... 

+sysn* [ .005818* (rand (ndp, 1) - . 5)  .02  90  9* (rand (ndp, 1) - . 5) 

.005818* (rand(ndp,l)-.5)  .01* (rand (ndp, 1) -. 5) ] ; 

%   PLOTS  FOR  VIEWING  ON  MONITOR  AND  STORED  IN  META  FILE 

%   ELVATOR  vs  TIME 

plot (t, (180/pi) *simdata (  :  ,  1)  )  ; 

xlabel ('Time  (seconds) ') ;ylabel ('Elevator  Input  (degrees)'); 

ans=[' title (''' ,typac, '   SIMULATED  INPUT'');']; 

eval (ans) ; pause 

%meta  A:\plots\deltae 

%  AOA  vs  TIME 

plot (t,  (180/pi) *simdata(: ,2) ) ; 

xlabel('Time  (seconds) ') ;ylabel ( 'AOA  Output  (degrees)'); 

ans=  [ 'title  (''' ,typac,  '   SIMULATED  DATA");']; 

eval (ans) ; pause 

Imeta  A:\plots\AOA 

%  PITCH  RATE  (Q)  vs  TIME 

plot (t,  (180/pi) *simdata ( : , 3) ) ; 

xlabel('Time  (seconds) ') ;ylabel (' Pitch  Rate,  Q,  Output  (deg/sec) ' ) ; 

ans=  ['title  ('"  ,typac,  '   SIMULATED  DATA");']; 

eval (ans) ; pause; 

Imeta  A:\plots\Q_ 

%  THETA  vs  TIME 

plot (t,  (180/pi) *simdata(: ,4) ) ; 

xlabel ('Time  (seconds) ' ) ;ylabel ('Theta  Output  (deg) ' ) ; 

ans=[' title  ('  '  '  ,typac,  '   SIMULATED  DATA");']; 

eval (ans) ; pause; 

%meta  A:\plots\theta 

%  NORMAL  ACC  vs  TIME 

plot (t, simdata ( : , 5) ) ; 

xlabel ('Time  (seconds) ') ;ylabel ('Normal  Acceration  Output  (G) ' ) ; 

ans=  [' title  ('"  ,typac,  '   SIMULATED  DATA");']; 

eval (ans) ; pause; 

%meta  A:\plots\G 

dispC  ') 

disp  (' '  )  ; 

dispC  ') 

disp ('NOTE  ====>   DATA  BEING  SAVED  TO  A   .mat  FILE') 
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dispC  FOR  MMLE  PROCESSING  ') 

if  sysn>0 

ans=['save  N', typac,'  typac  simdata  sysn  truth  iyy  gw  sref  cbar' ] ; 

eval (ans) 

disp('File  name  is  N  followed  by  the  type  aircraft .mat  N' )  , typac; 
else 

ans=['save  NN' , typac,'  typac  simdata  sysn  truth  iyy  gw  sref 
cbar'  ]  ; 

eval (ans) 

disp('File  name  is  NN  followed  by  the  type  aircraft. mat 
NN' ) , typac; 
end 

dispC  ') 
! dir/w 
dispC  ') 
dispC  ') 
dispC      NOW  RUN  npsmmlel.m  WITH  THE  CREATED  DATA  FILE.') 

disp  ( ' '  )  ; 

% END  LONGDAT .M 
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B 


LATDIR.M 


clear; 

%   MACRO  FILE  NAME 

% 


>=======  LATDIR.M  =======< 

Date:  3  Feb  92 


disp 
disp 
disp 
disp 
disp 
disp 


disp 
disp 
disp 
disp 
disp 
disp 
disp 
disp 


ans) 

'  ') 

'  MACRO  TO  GENERATE  SIMULATED  LATERAL-DIRECTIONAL  DATA' ) 

'    USING  AILERON  &  RUDDER  INPUTS  WITH  OR  W/O  NOISE') 

ans) 


') 


B747 


UAV 


') 
OTHER' ) 

■-=') 


GET  AIRCRAFT  TYPE  TO  BE  USED  - 

AIRCRAFT  TYPES  AVAILABLE 

'  ') 

'NAVION    A4D    F104A     JETSTAR 

'  ') 

'NOTE  ========>  PROGRAM  IS  CASE  SENSITIVE  < 

'  ') 

'SELECT   "OTHER"   TO  INPUT  DATA  FOR  USER  DEFINED  AIRCRAFT 

'  ') 

typac=input  ('TYPE  IN  DESIRED  A/C  FROM  THE  ABOVE  LIST.    ','s') 

%  DETERMINE  IF  NOISE  IS  WANTED 

dispC  ') 

sysn=input (' INPUT  A  1  TO  INCLUDE  NOISE  AND  A  ZERO  FOR  NO  NOISE 

ndp=301;% NUMBER  OF  DATA  POINTS  -  15  SEC 

dt=.05;% TIME  STEP  FOR  THE  DATA 

amp=.05/% AMPLITUDE  OF  INPUT  (RADIANS) 

period=l;% PERIOD  OF  DOUBLET  IN  SEC  =  PERIOD  +  1 

t=[0:ndp-l]  *dt;% TIME  VECTOR 

simdata=zeros (ndp,7) ;% SETUP  DATA  MATRIX  ALL  ZEROS 

% GENERATE  THE  INPUT  DOUBLET 

dslope= (4*amp*dt) /period; 

dl= (0 : -l*dslope : -l*amp) ; dla=-amp*ones (1:10) ; 

d2= (-l*amp+dslope : dslope : 0) ; 

d3= (dslope : dslope : amp) ;d3a=-l*dla; 

d4= (amp-dslope : -l*dslope : 0) ; 

% AILERON  AND  RUDDER  INPUTS 

simdata (:, 1)= [zeros (1, 60)  dl  dla  d2  d3  d3a  d4 

zeros (1, ndp- (period/dt) -81) ] ' ; 

simdata ( : ,2)=[-dl  -dla  -d2  -d3  -d3a  -d4 

zeros (1, ndp- (period/dt) -21) ] ' ; 

vtrue=input (' INPUT  AIRCRAFT  TRUE  VELOCITY     ft/sec 

altft=input (' INPUT  AIRCRAFT  ALTITUDE  ft 

oat =input (' INPUT  THE  OUTSIDE  AIR  TEMPERATURE  Deg  F 

if  (strcmpC  OTHER'  ,typac)>0)  ; 

typac=input (' INPUT  THE  A/C  TYPE  . . . .<  6  characters  ','s') 
sref=input (' INPUT  THE  AIRCRAFT  REFERENCE  AREA  ftA2  ') 
gw=input (' INPUT  THE  AIRCRAFT  GROSS  WEIGHT  lbs  ') 
ixx=input  ('AIRCRAFT  Ixx  MOMENT  OF  INERTIA    slug-ft/N2  '); 


') 


') 
') 
') 


'); 
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ixz=input ('AIRCRAFT  Ixz  MOMENT 
izz=input ('AIRCRAFT  Izz  MOMENT 
bbar=input ('AIRCRAFT  WING  SPAN 


OF  INERTIA 
OF  INERTIA 
LENGTH  "b" 


slug-ft 
slug-ft 
ft 


') 

') 


CYb=input (' INPUT  DERIVATIVE  CY_b 
Clb=input (' INPUT  DERIVATIVE  Cl_b 
CNb=input (' INPUT  DERIVATIVE  CN_b 
Clp=input (' INPUT  DERIVATIVE  Cl_p 
CNp=input (' INPUT  DERIVATIVE  CN_j> 
Clr=input (' INPUT  DERIVATIVE  Cl_r 
CNr=input (' INPUT  DERIVATIVE  CN  r 


CONTROL 
CONTROL 
CONTROL 
CONTROL 
CONTROL 


DERIVATIVE 
DERIVATIVE 
DERIVATIVE 
DERIVATIVE 
DERIVATIVE 


1/RAD 
1/RAD 
1/RAD 
1/RAD 
1/RAD 
1/RAD 
1/RAD 
Cl_da 
CN_da 
CY_dr 
Cl_dr 
CN  dr 


1/RAD 
1/RAD 
1/RAD 
1/RAD 
1/RAD 


MAT 
Cldr 


FILE 
CNdr] 


DATA'  ON  SELECTD  A/C 
CNda  CYdr  Cldr  CNdr]  , 


Clda=input (' INPUT 

CNda=input ( ' INPUT 

CYdr=input ( ' INPUT 

Cldr=input ( ' INPUT 

CNdr=input ( ' INPUT 

% 

ans=['save  ' , typac, ' .mat ; ' ] ; 

eval(ans);% SAVE  THE  NEW  A/C  DATA  IN  A 

truth=[CYb  Clb  CNb  Clp  CNp  Clr  CNr  Clda  CNda  CYdr 
else 

ans=['load  ' , typac, ' .mat ; ' ] ; 

eval(ans);% LOADS  THE  'TRUTH 

truth=[CYb  Clb  CNb  Clp  CNp  Clr  CNr  Clda 
end 

% CALCULATE  DENSITY,  DYNAMIC  PRESSURE,  AND  MASS 

rho=. 0023769*exp ( (-1*32 . 17*altft) / (1716* (oat+460) ) ) ; 
qbar= . 5*rho*vtrue*vtrue; 
mass=gw/32 . 17 ; 

% CALCULATE  CONSTANTS 

constl= (qbar*sref ) /mass; 

const2=qbar*sref *bbar; const2a=const2*bbar/ (2*vtrue) / 

const3=constl/32 . 17; 

% DRYDEN  TURBULENT  MODEL  CONSTANTS 

if  altft<1750 

lw=altft; 

else 

lw=1750;% SCALE  OF  TURBULENCE 

end 

sigw=05;% RMS  VALUE  OF  TURBULENCE  IN  FT/SEC 

lamda= (vtrue/lw) ;k= (3*sigw"2) *vtrue/ (pi*lw) ;beta=vtrue/ (sqrt (3) *lw) 

r 

%    SYSTEM  STATE  SPACE  MATRICES 
%    INERT IAL  MATRIX 

In=[l    0    0   0  0   0 

0  ixx  -ixz  0  0   0, 

0  -ixz  izz  0  0   0, 

0    0    0   1  0   0, 

0    0    0   0  10 
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0    0    0   0   0   1]; 
an=[ (constl*CYb/vtrue)   0   -1 
0; 

(const2*Clb)   (const2a*Clp) 
0; 

(const2*CNb)   (const2a*CNp) 
0; 

0  1 

0; 

0  0 


(32.17/vtrue)   ( const l*CYb/vt rue) 
(const2a*Clr)   0   (const2*Clb) 


(const2a*CNr) 
0 
0 
0 


0  0 

-2*lamda] ; 
% 
bn=[0       (constl*CYdr/vtrue) 

(const2*Clda)     (const2*Cldr) 

(const2*CNda)      (const2*CNdr) 

0  0 

0  0  (sqrt (k) /vtrue) 

0        0   (sqrt (k) /vtrue) * (beta-2*lamda) 
% 
a=inv (In) *an;b=inv (In) *bn; 


0  (const2*CNb) 

0         0 

0         0 

0  -(lamdaA2) 


0 
0 
0 
0 
0 
0] 


c=[l 

0 

0 

0 

const3*CYb 
% 
d=[0 

0 

0 

0 

0 

% 


0 
0 
0 
0 
0] 


0  0 

0  0 

0  0 

0  0 

const3*CYdr  0 


0 
0 
0 
0 
0] 


%  SIMDATA( 

%  SIMDATA( 

%  SIMDATA( 

%  SIMDATA( 


: , l)=AILERON  INPUT  (RAD)   SIMDATA ( : , 2) =RUDDER  INPUT  (RAD) 
:,3)=BETA  (RAD)      SIMDATA (:, 4) =ROLL  RATE  (RAD/SEC) 
:,5)=YAWRATE  (RAD/SEC)    SIMDATA (: ,  6) =ROLL  ANGLE  (RAD) 
: , 7 ) =LATERAL  ACC 

[phi,gam]=c2d(a,b, dt) ;  %  CONVERT  TO  DISCRETE 

% RAND  NUMBER  GENERATOR    MEAN=0      VARIANCE=1 

rand ( ' normal ' ) ; 

u= [simdata ( : , 1) , simdata ( : , 2) , sysn*rand (ndp, 1) , ones (ndp, 1) ] ; 

% — A  INPUTS 

[OUTY  OUTX]=dlsim(phi,gani,c,d,u)  ; 

simdata ( : , 3:7)=OUTY(: , 1:5) ; 

rand ( ' uniform' ) ; 

% ADD   THE   MEASUREMENT   NOISE 

simdata ( : , 3 : 7) =simdata ( : , 3:7) +  .... 
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sysn* [ .005818* (rand (ndp, l)-.5)  .02  909* (rand (ndp, 1) - . 5) 

.02909* (rand (ndp, l)-.5)  .005818* (rand (ndp, l)-.5) 

.01* (rand (ndp, 1)-. 5) ] ; 

%       PLOTS  FOR  VIEWING  ON  MONITOR  AND  STORED  IN  META  FILE 

%    AILERON  INPUT 

subplot (211) ;plot (t, (180/pi) *simdata(: , 1) ) ; 

xlabel (' Time  (seconds) ') ;ylabel ('Aileron  Input  (degrees)'); 

ans=[' title (' ' ' , typac, '   SIMULATED  INPUT' ');'] ;eval (ans) ; 

%   RUDDER  INPUT 

subplot (212) ;plot (t, (180/pi) *simdata ( : ,  2) ) ; 

xlabel ('Time  (seconds) ' ) ;ylabel ('Rudder  Input  (degrees) ' ) ; 

ans=['  title  ("'  , typac,  '   SIMULATED  DATA' ');'] ;eval (ans) ; 

pause/  %meta  A:\plots\Latinput 

%   SIDE  SLIP  (Beta) 

subplot (111) ;plot  (t,  (180/pi) *simdata(: ,  3)  )  ; 

xlabel ('Time  (seconds) ') /ylabel (' Side  Slip   ,  B,  Output  (deg)'); 

ans=[ ' title (''' , typac, '   SIMULATED  DATA' ');'] ;eval (ans) ; 

pause;  %meta  A:\plots\sslipout 

%   ROLL  RATE  (P) 

plot(t,  (180/pi) *simdata(: ,4) ) ; 

xlabel ('Time  (seconds) '); ylabel (' Roll  Rate   ,P,  Output  (deg/sec)'); 

ans=[' title  (' ' ' , typac, '   SIMULATED  DATA' ');'] ;eval (ans) ; 

pause;  %meta  A:\plots\rollrout 

%  YAW  RATE  (R) 

plot (t, (180/pi) *simdata ( : , 5) ) ; 

xlabel ('Time  (seconds) ') ; ylabel (' Yaw  Rate   ,R,  Output  (deg/sec)'); 

ans= ['title (''' , typac, '   SIMULATED  DATA' ');'] ;eval (ans) ; 

pause;  %meta  A:\plots\yawout 

%  ROLL  ANGLE  (phi) 

plot (t,  (180/pi) *simdata (  :  ,  6)  )  ; 

xlabel ('Time  (seconds) ') ; ylabel (' Roll  Angle   Output  (deg)'); 

ans=[' title (' ' ' , typac, '   SIMULATED  DATA' ');']; eval (ans) ; 

pause;  %meta  A:\plots\rollaout 

%  LATERAL  ACCERATION  (G) 

plot  (t, simdata ( : , 7) ) ; 

xlabel ('Time  (seconds) ') ; ylabel (' Lateral  Acceleration  (G) ' ) ; 

ans=[' title (' ' ' , typac, '   SIMULATED  DATA' ');'] ;eval (ans) ; 

pause;  %meta  A:\plots\gout 

clc; 

dispC '  )  ; 

dispC  ') 

disp('NOTE  ====>   DATA  BEING  SAVED  TO  A   .mat  FILE') 

dispC  FOR  MMLE  PROCESSING  ') 

if  sysn>0 

ans=['save  N', typac,'  typac  simdata  sysn  truth  ixx  izz  ixz  gw  sref 
bbar' ] ; 

eval (ans) 
else 

ans= [ ' save  NN', typac,'  typac  simdata  sysn  truth  ixx  izz  ixz  gw 
sref  bbar' ] ; 

eval (ans) 
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end 

dispC  ') 

! dir/w 

dispC   NOW  RUN  npsmmle2.m  WITH  THE  CREATED  DATA  FILE.') 

disp  (' '  ) 

% END  LATDIR . M 
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C.      SIMULATED  LONGITUDINAL  MMLE 

1.   NPSMMLE1.M 

%       MACRO    NAME  >======    NPSMMLEl.M   =======< 

%  Date:  3   Feb    92 

clear; 

.'erase   npsnunlel.log; 

diary   npsnunlel.log 

disp  (' '  )  / 

disp('  NPS  PARAMETER  IDENTIFICATION  MACRO  FILE  ') 

disp  ('FOR  SIMULATED  FLIGHT  DATA  USING  SIMPLFIED  SHORT  PERIOD  '.) 
dispC        LONGITUDINAL  STABILITY  AND  CONTROL  DERIVATIVES') 
dispC  ') 

dispC '  )  ; 

npsinitl  % RUN  INITIALIZATION  MACRO 

format  compact, clc 

load  npsinitl; 

global  sref  cbar  gw  iyy  vtrue  qbar  dt  all  ql  thl  anl  del; 

% INPUT  FLIGHT  TEST  DATA 

uydata=zeros (ndp, 6) ;%  ESTABLISH  DATA  MATRIX  FOR  MMLE.M 

% UYDATA ( 

% UYDATA  ( 

% UYDATA ( 

% UYDATA  ( 

% UYDATA  ( 

% COLUMN  NUMBER  ONE  IN  DATA  FILE 

coll= [ ' uydata ( : , 1) =simdata ( : , 1) ; ' ] ; eval 

% COLUMN  NUMBER  TWO  IN  DATA  FILE 

uydata ( :  ,  2) =ones (ndp, 1) ; 

% COLUMN  NUMBER  THREE  IN  DATA  FILE  ANGLE  OF  ATTACK 


,1)  - 

DELTA 

E 

(INPUT) 

,3)  = 

AOA 

,4)  = 

PITCH 

RATE  (Q) 

,5)  = 

THETA 

>(:,6) 

=  NORMAL 

ACC 

ELEVATOR  INPUT 

coll) , 

• 

UNITY 

INPUT 

col3= [ ' uydata ( : , 3) =simdata ( 

% COLUMN  NUMBER  FOUR 

col4= [ ' uydata ( : , 4 ) =simdata ( 

% COLUMN  NUMBER  FIVE 

col5= [ ' uydata ( : , 5) =simdata ( 


,2) ; ' ] ;eval (col3) ; 

IN  DATA  FILE  PITCH  RATE 

,3) ;'];eval(col4) ; 

IN  DATA  FILE  THETA 

,4);'];eval(col5); 

% COLUMN  NUMBER  SIX  IN  DATA  FILE  NORMAL  ACC 

col 6= [ ' uydata ( : ,  6) =simdata (:,5);'] ;eval (col 6) ; 

% INITIAL  CONDITIONS 

del=uydata (1,1) ;all=uydata (1,3) ;ql=uydata (1,4) ; 
thl=uydata (1,5) ;anl=uydata (1,6); 

% ADDITIONAL  INPUTS  TO  MMLE  FOLLOW 

p2snam='npsp2ssl' ;  %  MACRO  NAME  FOR  P2SS  FUNCTION 

pO=pref;  %-  INITIAL  PARAMETER  ESTIMATES  INPUT  DURING  NPSINIT1.M 

% —  CHECK  IF  THE  WEIGHTING  FUNCTION  IS  TO  BE  USED  FOR  INITIAL 

VALUES 

disp ('DO  YOU  WANT  TO  WEIGHT  THE  INITIAL  ESTIMATES?  INPUT  1=Y  0=NO 

') 

input  ( '  '); 

if  ans==l 
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disp ('Input  Weighting  Row  Vector  length  1x5') 

disp  ('Use  brackets-  ex.  [.111. 11]   &  lower  #  higher  weight') 

rmsO=input ('  ' ) ; 

end 

pidq=[l];% IDENTIFY  WHICH  PARAMETERS  ARE  TO  BE  IDENTIFIED 

pidm=[l:5];% IN  THE  QUADRATIC,  MARQUARDT,  AND  FINAL  STAGES. 

pidf=[l:5] ;%-  pidq,  pidm,  pidf  MUST  BE  VALID  EVEN  IF  NOT  USED 

opt=[0  5  10  10  .02  .005  .001  1]; 

if  sysn==0, opt (4) =0;end 

%-  DEFAULT  ITERATIONS  AND  CONVERGENCE  CRITERIA,  IF  NOISE  FREE 

OPT4=0 

gg0=eye  (4)  *  ( .  001)  ;% INNOVATIONS  COVARIANCE  MATRIX 

pert=le-4;%-PERTURBATION  USED  FOR  NUMERICAL  GRADIENT  CALCULATION 

linesearch=l;% USE  LINESEARCH  TO  HELP  PROC.  NOISE  CONVERGENCE 

!cls 

mmle  % CALL  MAIN  MMLE  MACRO  FROM  TOOL  BOX 

% PERFORM  FINAL  CALCULATIONS 

cla=pfin(l) ;cma=pfin (2) ;cmq=pf in (3) ;clde=pf in (4) ;cmde=pfin (5) ; 

dispC  ') 

deriv=[cla  cma  cmq  clde  cmde] ; 

dispC  MMLE    STABILITY  &  CONTROL  DERIVATIVES        ') 

dispC      CLA        CMA        CMQ        CLDE       CMDE') 

disp (deriv) 

dispC  INITIAL  INPUT  DERIVATIVES      ') 

disp (pref ) 

if  exist ('truth' ) ==1 ; 

dispC   TRUTH  DERIVATIVES  USED  TO  GENERATE  THE  DATA  '); 
disp (truth) 

end 

pause 

mleplotl 

diary  off 

%! print  npsmmlel.log; 

% END  NPSMMLE1.M 
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2.   NPSINIT1.M 

clear; 

%  MACRO  FILE  NAME  >=======  NPSINIT1.M  =======< 

%  Date:  3  Feb  92 

%  INITIAL  SETUP  MACRO  FOR  RUNNING  THE  PARAMETER 

%  ESTIMATION  TOOLBOX  IN  MATLAB  FOR  SIMPLIFIED  LONGITUDINAL 

%  SHORT  PERIOD  STABILITY  AND  CONTROL  DERIVATIVES 

%  THIS  MACRO  'GETS'  CONSTANTS  AND  DATA  FILE 

% LOAD  DATA  FILE 

dispC  DATA  FILE  MUST  CONTAIN  DATA  IN  COLUMN  MATRIX   ') 

disp('  MATRIX  NAME  SIMDATA:  N  DATA  PTS  X  5  COLUMNS  ') 

dispC  ELEVATOR/ALPHA/THETA/PITCH  RATE/NORMAL  ACC . ' ) 

dispC  ') 

dispC  DATA  FILE  NAME — MUST  EXIST  WITH  A   .mat   EXTENSION.  '); 

data=input (' ENTER  DATA  FILE  NAME  (  ==>  WITH  OUT  <==  .MAT 

EXTENSION)?   ','s'); 

if  exist ('dt' ) ==0, dt=input ('DELTA  T  BETWEEN  DATA  POINTS?  ' ) ; end 

ldc=['load  ', data, ' .mat; '] ; 

eval(ldc);% — EXECUTES  LOAD  COMMAND 

ndp=length (simdata) ; 

t=[0:ndp-l] *dt; 

% INPUT  REQUIRED  CONSTANTS IF  NOT  IN  DATA  FILE 

if  exist ('sref ' )==0, sref =input ('REFERENCE  AREA  (S)  IN  SQUARE  FEET? 

' ) ; end 

if  exist ('cbar' ) ==0, cbar=input ('MEAN  AERODYNAMIC  CHORD  IN  FEET? 

' ) ; end 

if  exist ('gw' ) ==0, gw=input ('AIRCRAFT  GROSS  WEIGHT  IN  POUNDS?  ' ) ; end 

if  exist (' iyy' )==0, iyy=input ('MOMENT  OF  INERTIA  (IYY)  IN  SLUG-FTA2? 

' ) ; end 

%sdc=['save   ',data,'   simdata  dt  sref  cbar  gw  iyy  ' ] ; eval (sdc) ; 

vtrue=input ('AIRSPEED  IN  FEET  PER  SECOND?  '); 

altft=input ('AIRCRAFT  ALTITUDE  IN  FEET?   '); 

oat=input ('OUTSIDE  AIR  TEMPERATURE?     '); 

% CALCULATE  CONSTANTS  DENSITY  AND  DYNAMIC  PRESSURE 

rho=.0023769*exp( (-1*32 . 174*altft) / (1716* (oat+460) ) ) ; 
qbar= . 5*rho*vtrue*vtrue; 

% INPUT  INITIAL  ESTIMATES  FOR  STABILITY 

% AND  CONTROL  DERIVATIVES 

pref (l)=input ('CL_ALPHA  ESTIMATE  (1/RAD) ?  '); 

pref (2) =input ('CM_ALPHA  ESTIMATE  (1/RAD)?  '); 

pref (3)=input ( ' CM_Q  ESTIMATE  (1/RAD)?  ') 

pref (4)=input ( ' CL_DE  ESTIMATE  (1/RAD)?  ' 

pref (5)=input ( ' CM_DE  ESTIMATE  (1/RAD)?  ' 

! erase  npsinitl .mat ; 

save  npsinitl 

% END  NPSINITL  M 
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3.   NPSP2SS1.M 


function  [a, phi, gam, c, d, q, xO, dt, rowinq,b] =npsp2ssl (p) 
MACRO  FILE  NAME     >=======  NPSP2SS1.M  =======< 

Date:  3  Feb  92 


MACRO  TO  EST.  FUNCTION  FOR  TRANSFORMING 
MODEL  PARAMETERS  INTO  STATE  SPACE  EQUATIONS 


P2SS  FUNCTION  FOR  NPSMMLE1.M 


P(D 
P(2) 
P(3) 
p(40 
P(5) 


CL 
CM" 

cm" 
cl" 
cm" 


ALPHA 
"ALPHA 

'Q 

DE 

DE 


STABILITY  AND  CONTROL 
PARAMETERS 


PERFORM  INITIAL  CALCULATIONS 

constl= (-l*qbar*sref ) / (cos (all) *gw*vtrue/32 . 17) ; 

const 2=qbar * sref*cbar/iyy; 

const3=const2*cbar/ (2*vtrue) / 

const4=32 . 17*cos (thl) / (vtrue*cos (all) ) ; 

const5= (qbar*sref ) /gw; 

const6=-l* (constl*p (1) *all+ql+constl*p (4) *del+const4) ; 

const7=-l* (const2*p (2) *all+const3*p (3) *ql+const2*p (5) *del) ; 

const8= (-1* (const5*p (1) *all+const5*p (4) *del) ) +anl; 

% STABILITY  DERIVATIVES 

a= [constl*p  (1)  1  0; 

const2*p(2)  const3*p(3)       0; 

0  10]; 

% CONTROL  DERIVATIVES 

b= [constl*p (4)      (const4+const6) / 

const2*p(5)  const7; 

0  (-l*ql)]; 

% 

c=[l  0 

0  1 

0  0 

const5*p(l)  0 

% 

d= [ 0  0 , 


-MEASUREMENT 

0, 

0, 

1 
0] 
FEED 


MATRIX 


THROUGH  MATRIX 


0 
0 
const5*p (4) 


0 

0 

const8] 


%■ 

q=eye (a) *le-4 ; % — 

% 


rowinq=[0    0    0    0 


Q    IS 

-  WITH 

-  ROWS 
SAME 

0]/ 


STATE  NOISE  COVARIANCE 

THE  SAME  SIZE  AS  a 

Q*Q'  POS.  DEFINITE 

IN  Q  IN  WHICH  PARAMETERS  OCCUR,  A  VECTOR 

DIMENSION  AS  p 
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% INITIAL  STATE  VECTOR 

xO=[all  ql  thl] ; 

% DISCRETIZE 

%  *****NEED  TO  EDIT  dt  (below)  FOR  THE  DELTA  T  OF  THE  DATA  ***** 

dt=.05; 

[phi, gam] =c2d (a,b, dt) ; 

% END  NPSP2SS1.M 
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4.   MLEPL0T1.M 

%  MACRO  FILE  NAME         >=======  MLEPLOT1.M  =======< 

%  Date:  3  Feb  92 

%  MACRO  TO  PLOT  DATA  FROM  NPSMMLE1 

constl= (-l*qbar*sref ) / (cos (all) *gw*vtrue/32 . 17) ; 

const 2=qbar * s re f*cbar/iyy; 

const3=const2*cbar/ (2*vtrue) ; 

const4=32 . 17*cos (thl) / (vtrue*cos (all) ) ; 

const5= (qbar*sref ) /gw; 

const6=-l* (constl*pf in (1) *all+ql+constl*pf in (4) *del+const4) ; 

const7=-l* (const2*pf in (2) *all+const3*pf in (3) *ql+const2*pf in (5) *del) 


const8=  (-1*  (const5*pfin (1) *all+const5*pf in (4) *del) ) +anl; 

% STABILITY  DERIVATIVES 

a=[constl*pfin  (1)  1  0, 

const2*pfin(2)  const3*pf in (3)  0 

0  10] 

% CONTROL  DERIVATIVES 

b= [constl*pf in  (4)      (const4+const6) ; 

const2*pf in (5)  const7; 

0  (-l*ql)]; 

% MEASUREMENT  MATRIX 

c=[l  0  0, 

0  10 

0  0  1 

const5*pfin  (1)  0  0] 

% FEED  THROUGH  MATRIX 

d=  [  0  0 ; 

0  0; 

0  0; 

const5*pf in (4)     const8]; 
rtdc=(180/pi) ; 

% OUT2  =  OUTPUT  VECTOR OUT3  =  STATE  VECTOR 

[OUT2,OUT3]=lsim(a,b, c,d,uydata( : , 1 :2) , t,x0) ; 

if  exist ('truth' ) >0; 

const 6=-l* (constl*truth (1) *all+ql+constl*truth (4) *del+const4) ; 

const7=-l* (const2*truth (2) *all+const3*truth (3) *ql+const2*truth (5) *d 

el)  ; 

const8= (-1* (const5*truth(l) *all+const5*truth (4) *del) ) +anl; 


a=[constl*truth(l) 
const2*truth(2) 
0 


b=[constl*truth(4) 

const2*truth(5) 

0 
% 

c=[l 

0 


1  0; 

const3*truth(3)       0; 
1  0]; 

CONTROL  DERIVATIVES 

(const4+const6) ; 
const7; 
(-l*ql)]; 
MEASUREMENT  MATRIX 

0  0; 

1  0; 
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0  0  1; 

const5*truth(l)  0  0] ; 

% FEED  THROUGH  MATRIX 

d=[0  0; 

0  0; 

0  0; 

const5*truth (4)     const8]; 
[TRU2,TRU3]=lsim(a,b,c,d,uydata(: ,1:2) ,t,x0) ; 
end 

% PLOTS  TO  MONITOR  AND  STORED  IN  META  FILE 

if  exist ('typac' ) ==0, typac=input (' INPUT  THE  AIRCRAFT  TYPE  ? 

' , ' s ' ) ; end 

.'erase  a:\plots\outputg.met; 

! erase  a:\plots\outputth.met; 

.'erase  a:\plots\outputde.met; 

.'erase  a:\plots\outaoa.met; 

.'erase  a:\plots\outputq.met; 

% ELEVATOR  VS  TIME 

hold  of f ;plot (t, rtdc*uydata ( : , 1) ,  '  -r'  ) ; 

xlabel('Time  (seconds) ') ;ylabel (' Elevator  Input  (degrees)'); 

ans=[' title  (" ' , typac, '  ELEVATOR  INPUT  VS  TIME  '');']; 

eval (ans) ; pause 

meta  A:\plots\OUTPUTde 

% AOA  (OBSERVED  AND  ESTIMATED)  VS  TIME 

plot (t, rtdc*uydata ( : , 3) , ' *r' ) ;hold  on; 

xlabel('Time  (seconds) ') ;ylabel ('AOA  (degrees)'); 

ans=[' title (' ' ' , typac, '   ESTIMATED  AND  MEASURED  AOA  RESPONSE'');'] 

eval (ans) 

text ( . 6, . 85, ' *   Measured  Data  Points  ' , ' sc' ) /pause; 

plot  (t , rtdc*OUT2 ( : , 1 ) , ' og'  )  ; 

text ( . 6, . 80, ' o   Estimated  Response  ' , ' sc' ) ;pause; 

if  exist ('truth' ) >0; 

plot (t, rtdc*TRU2 ( : , 1) , ' -b' ) ; 

text ( . 6, . 75, ' -   True  Response  ' , ' sc' ) ; pause; 
end 
pause 
meta  A:\plots\outAOA 

% Q  (OBSERVED  AND  ESTIMATED)  VS  TIME 

hold  off;plot (t,rtdc*uydata(: ,4) , ' *r' ) ;hold  on; 

xlabel('Time  (seconds) ') ;ylabel (' Pitch  Rate,  Q,  (deg/sec) ' ) ; 

ans=[' title (' ' ' , typac, '   ESTIMATED  AND  MEASURED  q  RESPONSE'');']; 

eval (ans) 

text  (  . 6,  . 85, ' *   Measured  Data  Points  ' , ' sc' ) ;pause; 

plot (t, rtdc*OUT2 ( : , 2) , ' og' ) ; 

text  (  . 6,  . 80, ' o   Estimated  Response  ' , ' sc' ) ;pause; 

if  exist ('truth' ) >0; 

plot (t, rtdc*TRU2 ( : , 2) , ' -b' ) ; 

text ( . 6, . 75, ' -   True  Response  ' , ' sc' ) ; pause; 
end 
pause 
meta  A:\plots\OUTPUTQ 
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% THEATA  (OBSERVED  AND  ESTIMATED)  VS  TIME 

hold  off /plot (t, rtdc*uydata(: ,5) , '  *r' ) ;hold  on; 

xlabel (' Time  (seconds) ') ;ylabel ('Pitch  Angle,  Theta,  (deg) ' ) ; 

ans=[' title  (' ' ' ,typac, '   ESTIMATED  AND  MEASURED  THETA 

RESPONSE' ');']; 

eval (ans) 

text ( . 6,  . 85,  '  *   Measured  Data  Points  ','sc') ;pause; 

plot (t, rtdc*0UT2 ( : , 3) , ' og'  )  ; 

text ( . 6, . 80, ' o   Estimated  Response  ','sc') ;pause; 

if  exist ('truth' ) >0; 

plot (t,rtdc*TRU2 (: ,3) , '-to' ) ; 

text ( . 6, . 75, ' -   True  Response  ',' sc' ) /pause; 
end 
pause 
meta  A:\plots\OUTPUTTH 

% ACCELERATION  (OBSERVED  AND  ESTIMATED)  VS  TIME 

hold  of f ;plot (t, uydata ( : , 6) , ' *r' ) ;hold  on; 

xlabel ('Time  (seconds) ') ;ylabel ('Acceleration,  G  '); 

ans=[ ' title ("' ,typac, '   ESTIMATED  AND  MEASURED  G  RESPONSE'');']; 

eval (ans) 

text ( . 6, . 85, ' *   Measured  Data  Points  ' , ' sc' ) ;pause; 

plot(t,OUT2(:,4) , ' og' ) ; 

text  ( . 6,  . 80, ' o   Estimated  Response  ' , ' sc' ) ;pause; 

if  exist ('truth' ) >0; 

plot (t,TRU2 (: ,4) , ' -b' ) ; 

text ( . 6, . 75, ' -   True  Response  ' , ' sc' ) ;pause; 
end 
pause 

meta  A:\plots\OUTPUTG 
hold  off; 
% END  MLEPLOT1  .M 
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D.      SIMULATED  LATERAL-DIRECTIONAL  MMLE 
1.   NPSMMLE2.M 


MACRO  NAME 


% 

% 

clear; 

.'erase  npsmmle2.log; 

diary  npsmmle2 . log 


>======  NPSMMLE2.M  =======< 

Date:  3   Feb   92 


STABILITY  PARAMETER  IDENTIFICATION  MACRO  FILE') 

FOR  LATERAL-DIRECTIONAL  ' ) 

STABILITY  AND  CONTROL  DERIVATIVES  ' ) 


disp (ans) 
disp ('NPS 
disp  (' 
disp  (' 
dispC  ') 
disp (ans) 

npsinit2  % RUN  INITIALIZATION  MACRO 

format  compact , clc; 

load  npsinit2; 

global  sref  bbar  gw  ixx  ixz  izz  vtrue  qbar  dt  betal  rolll  yawl 

ranglel  ayl  dal  drl; 

% 

uydata=zeros (ndp, 8) ; %  


INPUT  FLIGHT  TEST  DATA 
ESTABLISH  DATA  MATRIX  FOR  MMLE . M 


%  UYDATA( 

(P) 

%  UYDATA( 
%  UYDATA( 
(phi) 

%  UYDATA( 
(ay) 


1)  =  DELTA  Aileron  (INPUT) 


UYDATA ( : , 5 )  =  ROLL  RATE 


,2) 

,3) 


=  DELTA 
=  UNITY 


Rudder 
INPUT 


(INPUT) 


UYDATA ( 
UYDATA ( 


,6) 
,7) 


=  YAW  RATE  (r) 
=  ROLL  ANGLE 


4)  =  BETA  (SIDE  SLIP) 


UYDATA ( : , 8 )  =  LATERAL  G 


coll  = 

% 

col2  = 


uydata(: , 3)=ones (ndp, 1) ;%  UNITY  INPUT 


NUMBER  ONE  IN 
: , 1) =simdata ( 


COLUMN 
' uydata ( 

COLUMN  NUMBER  TWO  IN 
' uydata {: ,2) =simdata ( 


DATA  FILE  AILERON  INPUT 
,1) ;'];eval(coll) ; 
DATA  FILE  RUDDER  INPUT 
,2);'];eval(col2); 


DATA  FILE  BETA  (SIDE  SLIP) 
,3);'];eval(col4); 


(P) 


% COLUMN  NUMBER  FOUR  IN 

col4= [ ' uydata ( : , 4 ) =simdata ( : 

% COLUMN  NUMBER  FIVE  IN  DATA  FILE  ROLL  RATE 

col5= [ 'uydata ( : , 5) =simdata ( : , 4) ; ' ] ;eval (col5) ; 

% COLUMN  NUMBER  SIX  IN  DATA  FILE  YAW  RATE  (r) 

col6= [ ' uydata (:,6) =simdata ( : , 5) ; ' ] ; eval (col6) ; 

% COLUMN  NUMBER  SEVEN  IN  DATA  FILE  ROLL  ANGLE 

col7= [ ' uydata ( : , 1) =simdata ( : , 6) ; ' ] ; eval (col7) ; 

% COLUMN  NUMBER  EIGHT  IN  DATA  FILE  LATERAL  G 

col8= [ ' uydata ( : , 8) =simdata ( : , 7) ; ' ] ;eval (col8) ; 

% INITIAL  CONDITIONS  FOR  da,  dr,  BETA,  p,  r,  phi, 

betal=uydata (1,4) ;rolll=uydata (1,5) ;yawl=uydata (1,6); 
ranglel=uydata (1,7) ;ayl=uydata (1,8)  ; 
dal=uydata (1,1) ;drl=uydata (1,2) ; 


(phi) 
(ay) 


ay 
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% ADDITIONAL  INPUTS  TO  MMLE  FOLLOW 

p2snam='npsp2ss2' ;% :-  MACRO  NAME  FOR  P2SS  FUNCTION 

pO=pref;% INITIAL  PARAMETER  ESTIMATES 

%rmsO= IF  USED  IT  IS  THE  WEIGHTING  FUNCTION  FOR  INITIAL  VALUES 

pidq=[l]/% IDENTIFY  WHICH  PARAMETERS  ARE  TO  BE  IDENTIFIED 

pidm=[l:12] ;% IN  THE  QUADRATIC,  MARQUARDT,  AND  FINAL 

STAGES  . 

pidf=[l :12] ;% pidq,  pidm,  pidf  MUST  BE  VALID  EVEN  IF  NOT  USED 

opt=[0  5  5  10  .02  .05  .001  1];%-  DEFAULT  ITERATIONS  AND  CONVERGENCE 
%  CRITERIA,  IF  NOISE  FREE  OPT4=0 

if  sysn==0, opt (4) =0;end 

gg0=eye  (5)  *  ( .  01)  ;  % INNOVATIONS  COVARIANCE  MATRIX 

pert=le-4;% —  PERTURBATION  USED  FOR  NUMERICAL  GRADIENT  CALCULATION 

linesearch=l;% USE  LINESEARCH  TO  HELP  PROC .  NOISE  CONVERGENCE 

mmle% CALL  MAIN  MMLE  MACRO  FROM  TOOL  BOX 

% PERFORM  FINAL  CALCULATIONS 

CYb=pfin (1) ;Clb=pfin (2) ;CNb=pfin(3) ;Clp=pfin(4) ; 
CNp=pfin(5) ;Clr=pfin(6) ;CNr=pfin(7) ;Clda=pf in (8) ; 
CNda=pfin(9) ;CYdr=pf in (10) ;Cldr=pf in (11) ;CNdr=pf in (12) ; 
deriv=[CYb        Clb         CNb        Clp        CNp        Clr; 

CNr        Clda        CNda       CYdr       Cldr       CNdr]  ; 
disp('  MMLE    STABILITY  &  CONTROL  DERIVATIVES        ') 

disp('CY_b        Cl_b        CN_b        Cl_p        CN_p      Cl_r  ') 
disp (deriv  (1,  :  ) ) 

disp('CN_r        Clda        CNda        CYdr        Cldr       CNdr') 
disp (deriv (2, : ) ) 

disp('  INITIAL  INPUT  DERIVATIVES      ') 

disp (pref (1:6)) ; 
disp (pref (7:12)  )  ; 
if  exist ('truth' ) ==1 / 

dispC         "TRUTH  DERIVATIVES"  USED  TO  GENERATE  DATA  ') 

disp (truth (1:6) ) ; 

disp (truth (7: 12)  )  ; 
end 

pause;mleplot2 
diary  off 

%  .'print  npsmmle2.log; 
% END  NPSMMLE2.M 
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NPSINIT2  M 


clear; 

%  MACRO  FILE  NAME  >=======  NPSINIT2.M  =======< 

%  Date:  3   Feb  92 

%  INITIAL  SETUP  MACRO  FOR  RUNNING  THE  PARAMETER 

%  ESTIMATION  TOOLBOX  IN  MATLAB  FOR  LATERAL-DIRECTIONAL 

%  STABILITY  AND  CONTROL  DERIVATIVES 

%  THIS  MACRO  'GETS'  CONSTANTS  AND  DATA  FILE 

disp('  ') 

dispC  DATA  FILE  NAME — MUST  EXIST  WITH  A   .mat   EXTENSION.   '); 

data=input ('ENTER  DATA  FILE  NAME  (  ==>  WITH  OUT  <==  .MAT 

EXTENSION) ?   '  ,  '  s'  )  ; 

ldc= [' load  ' , data, ' .mat; ' ] ;eval (ldc) ; 

ndp= length (simdata) ; 

if  exist ('dt' ) ==0, dt=input (' INPUT  dt  BETWEEN  DATA  POINTS.   ' ) ; end 

t=[0:ndp-l]*dt; 

% INPUT  REQUIRED  CONSTANTS 

if  exist (' sref ) ==0, sref=input ('REFERENCE  AREA  (S)  IN  SQUARE  FEET? 

'  )  ;  end 

if  exist ('bbar' ) ==0,bbar=input ('WINGSPAN  IN  FEET?  ');end 

if  exist ('gw' )==0,gw=input ('AIRCRAFT  GROSS  WEIGHT  IN  POUNDS?  ');end 

if  exist  ('  ixx'  )==0,  ixx=input  ('MOMENT  OF  INERTIA  (Ixx)  IN  SLUG-FT/N2? 

'  )  ;  end 

if  exist ('ixz' )==0,ixz=input ('MOMENT  OF  INERTIA  (Ixz)  IN  SLUG-FTA2? 

'  )  ;  end 

if  exist (' izz' )==0, izz=input ('MOMENT  OF  INERTIA  (Izz)  IN  SLUG-FTA2? 

' ) ; end 

%sdc=['save   '  , data, '   simdata  dt  sref  bbar  gw  ixx  ixz  izz 

' ] ; eval (sdc) ; 

vtrue=input ('AIRSPEED  IN  FEET  PER  SECOND?  '); 

altft=input ('AIRCRAFT  ALTITUDE  IN  FEET?   '); 

oat=input ('OUTSIDE  AIR  TEMPERATURE?     '); 

rho=.0023769*exp( (-1*32 . 174*altft) / (1716* (oat+460) ) ) ; 

qbar= . 5*rho*vtrue*vtrue; 
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% INPUT  INITIAL  ESTIMATES  FOR  STABILITY 

% AND  CONTROL  DERIVATIVES  TO  START  MMLE  PROGRAM 

pref (l)=-.6;%  input (' CY_beta  FROM  WIND  TUNNEL  (1/RAD) ? 

pref (2)=-. 15;%  input (' Cl_beta  FROM  WIND  TUNNEL 

pref (3)=. 20;%  input (' CN_beta  FROM  WIND  TUNNEL 

pref (4)=-. 35;%  input ('Cl_p     FROM  WIND  TUNNEL 

pref (5)=-. 05;%  input('CN_p     FROM  WIND  TUNNEL 

pref (6)=. 15;%  input ('Cl_r     FROM  WIND  TUNNEL 

pref (7)=-. 2;%  input ('CN_r     FROM  WIND  TUNNEL 

pref (8)=. 05;%  input (' Cl_da    FROM  WIND  TUNNEL 

pref (9)=-. 001;%  input('CN_da    FROM  WIND  TUNNEL  (1/RAD)?  '); 

pref (10)=.175;%input ('CY_dr    FROM  WIND  TUNNEL  (1/RAD)?  '); 

pref (ll)=.02;%input ( ' Cl_dr    FROM  WIND  TUNNEL  (1/RAD)?  ' ) ; 

pref (12)=-.075;%input ( ' CN_dr    FROM  WIND  TUNNEL  (1/RAD)?  '); 

.'erase  npsinit2  .mat; 

save  npsinit2 

% END  NPSINIT2.M 


(1/RAD) ? 

/   r 

') 

(1/RAD) ? 

'); 

(1/RAD) ? 

') 

(1/RAD) ? 

') 

(1/RAD) ? 

'); 

(1/RAD) ? 

'  ); 

(1/RAD) ? 

' ); 
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NPSP2SS2.M 


[a, phi, gam, c, d, q, xO, dt , rowinq,b] =npsp2ss2 (p) 
FILE  NAME     >=======  NPSP2SS2.M  =======< 

Date:  3   Feb  92 


function 
%   MACRO 

% 

% 

% 
% 


MACRO  TO  ESTABLISH  FUNCTION  FOR  TRANSFORMING 
MODEL  PARAMETERS  INTO  STATE  SPACE  EQUATIONS 


%  P2SS 

a. 

FUNCTION  FOR  NPSMMLE2.M 

%  P(D 

=  CY  beta   |                      |  p(7)  = 

CN  r 

%  p(2) 

=  CI  beta   | STABILITY  AND  CONTROL  |  p(8)  = 

CI  da 

%  p(3) 

=  CN  beta   |                       |  p(9)  = 

CN  da 

%  P(4) 

=  Cl_p       |       PARAMETERS         |  p(10)  = 

CY  dr 

%  p(5) 

=  CN_j>      |                         |  p(ll)  = 

CI  dr 

%  p(6) 

=  CI  r      |                        |  p(12)  = 

CI  dr 

% 

PERFORM  INITIAL  CALCULATIONS 

constl= (qbar*sref ) / (gw/32 . 17)  ; 

const2=qbar*sref *bbar; const2a=const2*bbar/ (2*vtrue) ; 

const3= (qbar*sref ) /gw; 

% INERTIAL  MATRIX 

In=[l   0      0   0 

0   ixx  -ixz  0 

0  -ixz   izz  0 

0    0     0   1]; 

% PLANT 

an= [constl*p (1) /vtrue  0 

const2*p(2)       const2a*p(4) 

const2*p(3)       const2a*p(5) 

0  1 

% 

a=inv (In) *an; 
% 
bn=[0        constl*p (10) /vtrue 


-1 


(32.17/vtrue) 


const2a*p (6) 

const2a*p (7) 

0 


0 
0, 
0] 


const2*p(8) 
const2*p  (9) 
0 

% 

b=inv (In) *bn; 


const2*p(ll) 
const2*p  (12) 
0 


0; 
0; 
0; 

0] 


c=[l 

0 

0 

0, 

0 

1 

0 

0, 

0 

0 

1 

0, 

0 

0 

0 

1, 

const3*p  (1) 

0 

0 

0], 

% 

d=[0          0 

betal; 

0           0 

rolll; 

0           0 

\ 

f  awl ; 

81 


0  0  r angle 1; 

0  const3*p(10)       ayl]; 

% 

% STATE  NOISE  COVARIANCE 

q=eye(a)  *le-4;% Q  IS  THE  SAME  SIZE  AS  a 

%  WITH  Q*Q'  POS.  DEFINITE! 

% ROWS  IN  Q  IN  WHICH  PARAMETERS  OCCUR,  A  VECTOR 

%  SAME  DIMENSION  AS  p 

rowinq=0*p; 

% INITIAL  STATE  VECTOR 

xO=[betal  rolll  yawl  ranglel]; 

% DISCRETIZE 

dt=.05; 

[phi, gam] =c2d (a,b, dt) ; 

% END  NPSP2SS2.M 
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MLEPL0T2 . M 


%  MACRO  FILE  NAME    >=======  MLEPLOT2.M  =======< 

%  Date:  3   Feb   92 

% MACRO  TO  PLOT  DATA  FROM  NPSMMLE2 

% GENERATE  THE  PREDICTED  DATA 

% CALCULATE  STABILITY  AND  CONTROL  MATRIX 

mass=gw/32 . 17/ rtdc= (180/pi) ; constl= (qbar*sref ) /mass; 
const2=qbar*sref *bbar; const2a=const2*bbar/ (2*vtrue) ; 
const 3=qbar*sref/gw; 


In=[l    0    0    0; 

0   ixx  -ixz  0; 

0  -ixz   izz  0; 

0    0    0    1] 

an= [constl*pf in (1) /vt 

rue        0              -1       32 . 

1 7 /vtrue; 

const2*pfin(2) 

const2a*pf in (4)   const2a*pf in (6) 

0; 

const2*pfin (3) 

const2a*pf in (5)   const2a*pf in (7) 

0; 

0 

1                0 

0] 

a=inv (In) *an; 

% 

bn=[0 

constl*pf in (10) /vtrue           0, 

const2*pfin (8) 

const2*pfin (11)           0, 

const2*pfin (9) 

const2*pfin(12)           0, 

0 

0          0]  , 

b=inv (In) *bn; 

% 
c=[l 

0 

0       0; 

0 

1 

0       0; 

0 

0 

1       0; 

0 

0 

0       1; 

const3*pf in (1) 
% 
d=[0          0 

0 

0      0]; 

betal; 

0          0 

rolll; 

0          0 

yawl ; 

0           0 

ranglel ; 

const3*pf in (10) 


ayl]; 


%  OUT2  =  OUTPUTS  OUT3  =  STATE  VECTOR 

[OUT2,OUT3]=lsim(a,b, c,d,uydata ( : , 1 : 3) , t , xO) ; 

if  exist ( ' truth' ) >0 

an=[constl*truth(l) /vtrue      0  -1       32.17/vtrue; 

const2*truth(2)       const2a*truth (4)   const2a*truth (6) 
0; 

const2*truth(3)       const2a*truth (5)   const2a*truth (7) 
0; 

0  1  0  0  ]  ; 

a=inv (In) *an; 
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constl*truth (10) /vtrue 
const2*truth(ll) 
const2*truth(12) 
0 


0; 
0/ 
0; 
0]; 


0 
0 
0 

1 

0]; 


0 
0 
0 
0 
const3*truth(10) 


0 
0 
0 
0 
0] 


% 

bn=[0 

const2*truth(8) 
const2*truth(9) 
0 
b=inv(In) *bn; 
% 
c=[l 

0 

0 

0 

const3*truth (1) 
% 
d=[0 

0 

0 

0 

0 

% 

%TRU2  =  OUTPUT   TRU3  =  STATE  VECTOR 
[TRU2,TRU3]=lsim(a,b,c,d,uydata( : , 1:3) ,t) ; 

end 

%     PLOTS  FOR  VIEWING  ON  MONITOR  AND  STORE  TO  META  FILE 

! erase  a : \plots\* .met 

subplot (211) ;plot (t,rtdc*uydata(:  ,1)); 

xlabel ( ' Time  (seconds) ') ;ylabel ( 'Aileron  Input 

ans=  ['  title  ('"  , typac,  '  AILERON  INPUT  VS  TIME 

subplot (212) ;plot (t , rtdc*uydata ( :  ,  2)  )  / 

xlabel ('Time  (seconds) ') ;ylabel (' Rudder  Input 

ans= [ 'title (''' , typac, '  RUDDER  INPUT  VS  TIME  ' 

pause;  %meta  A:\plots\INPUTDAR 

%   Beta  vs  Time 

subplot  (111)  ;plot  (t,  rtdc*uydata  (:,4)/*r')  ;hold  on; 

xlabel ('Time  (seconds) ') ;ylabel (' Beta  (degrees)'); 

ans= [ 'title (''' ,typac, '   ESTIMATED  AND  MEASURED  Beta 

RESPONSE' ');']; 

eval (ans) 

text ( . 6, . 85, ' *   Measured  Data  Points 

plot (t,rtdc*OUT2 (: , 1) , ' og' ) ; 


(degrees) ' ) ; 
'');'] ;eval (ans) 

(degrees) ' ) ; 
' ) ; ' ] ; eval (ans) 


text  ( . 6,  . 80, ' o   Estimated  Response 


' , ' sc' ) ; pause; 
' sc' ) ; pause; 


if  exist ('truth' ) >0; 

plot (t,rtdc*TRU2 (: , 1) , '-b' ) ; 

text ( . 6, . 75, ' -   "True  Response"   ','sc'); 

end/pause;  %meta  A:\plots\outbeta 
%  Roll  rate  vs  Time 

hold  of f ;plot (t, rtdc*uydata ( : , 5) , ' *r' ) ;hold  on; 
xlabel ('Time  (seconds) ') ;ylabel (' Roll  Rate,  p,  (deg/sec) ' ) ; 
ans=[' title (' ' ' ,typac, '   ESTIMATED  AND  MEASURED  p  RESPONSE'') 
eval (ans) 

text ( . 6, . 85, ' *   Measured  Data  Points  ' , ' sc' ) ;pause; 
plot (t,rtdc*OUT2 (:,2) , ' og' ) ; 


'] 
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text ( . 6, . 80, ' o   Estimated  Response  ' , ' sc' ) ;pause; 

if  exist ('truth' ) >0; 

plot (t, rtdc*TRU2 ( : , 2) , ' -b' ) ; 

text  ( .  6,  .  75,  '  -   "True  Response"  '/sc'); 

end/pause; %meta  A:\plots\OUTP 
%   Yaw  Rate  vs  Time 

hold  off ;plot (t, rtdc*uydata( : , 6) , ' *r' ) ;hold  on; 
xlabel('Time  (seconds) ') ;ylabel (' Yaw  Rate,  r,  (deg/sec) ' ) ; 
ans=[' title (' ' ' ,typac, '   ESTIMATED  AND  MEASURED  r  RESPONSE'');']; 
eval (ans) 

text ( . 6,  . 85,  '  *   Measured  Data  Points  ' , ' sc' ) ;pause; 
plot  (t, rtdc*OUT2 ( : , 3) , ' og' ) ; 
text ( . 6,  . 80,  ' o   Estimated  Response  ' , ' sc' ) ; pause; 

if  exist ('truth' ) >0; 

plot (t, rtdc*TRU2 ( :  ,  3)  ,  '  -b'  )  ; 

text ( . 6, . 75, '-   "True  Response"  ','sc'); 

end;pause; %meta  A:\plots\OUTr 
%  Bank  Angle  vs  Time 

hold  off ;plot (t,rtdc*uydata ( : ,7) , ' *r' ) ;hold  on; 
xlabel('Time  (seconds) ') ;ylabel (' Bank  Angle,  phi,  (deg) ' ) ; 
ans= [ 'title ('"  ,typac, '   ESTIMATED  AND  MEASURED  phi  RESPONSE'');']; 
eval (ans) 

text ( . 6,  . 85,  '  *   Measured  Data  Points  ' , ' sc' ) ;pause; 
plot (t, rtdc*0UT2 < :  ,  4) , ' og' ) ; 
text ( . 6, . 80, ' o  Estimated  Response  ', 'sc') /pause; 

if  exist ('truth' ) >0; 

plot (t,rtdc*TRU2 ( : , 4) , '-b' ) ; 

text ( . 6,  .75,  '  -   "True  Response"  ','sc'); 

end;pause;  %meta  A:\plots\OUTphi 
%  Lateral  G  vs  Time 

hold  off ;plot (t,uydata(: , 8) , ' *r' ) ;hold  on; 
xlabel('Time  (seconds) ') ;ylabel (' Lateral  G,  ay,  (G) ' ) ; 
ans=[' title (' ' ' ,typac, '   ESTIMATED  AND  MEASURED  Lateral  G 
RESPONSE' ');']; 
eval (ans) 

text ( . 6,  . 85,  '  *   Measured  Data  Points  ' , ' sc' ) ;pause; 
plot (t,OUT2 (: ,5) ,  '  og'  )  ; 
text ( . 6,  . 80,  '  o   Estimated  Response  ' , ' sc' ) ;pause; 

if  exist ('truth' ) >0; 

plot (t,TRU2 (: ,5) , '-b' ) ; 

text  (  . 6,  . 75,  '  -   "True  Response"  ','sc'); 

end;pause;  %meta  A:\plots\OUTlatg 
hold  off; 
% END  MLEPLOT2 .M 
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E.      ACTUAL  LONGITUDINAL  MMLE 
1.   NPSMMLE3.M 


%   MACRO  NAME         >======  NPSMMLE3.M  =======< 

%  Date:  31  Jan  91 

clear; 

!  erase  npsnunle3.log; 

diary  npsmmle3.log 

disp  (' '  )  ; 

dispC       NPS  PARAMETER  IDENTIFICATION  MACRO  FILE  ') 

disp('    FOR  ACTUAL  FLIGHT  TESTS  USING  SIMPLFIED  SHORT  PERIOD    ') 

dispC        LONGITUDINAL  STABILITY  AND  CONTROL  DERIVATIVES') 

dispC  ') 

disp  (' '  )  ; 

npsinit3  % RUN  INITIALIZATION 

MACRO 

format  compact, clc 

load  npsinit3; 

global  sref  cbar  gw  iyy  vtrue  qbar  dt  all  ql  thl  anl  del; 

% INPUT  FLIGHT  TEST  DATA 

uydata=zeros (ndp, 6) ;%  ESTABLISH  DATA  MATRIX  FOR  MMLE.M 

% UYDATA ( 

% UYDATA ( 

% UYDATA  ( 

% UYDATA ( 

% UYDATA  ( 

% COLUMN  NUMBER  ONE  IN  DATA  FILE  ELEVATOR  INPUT 

coll= [ ' uydata ( : , 1) =simdata ( : , 1) ; ' ] ; eval (coll) ; 

% COLUMN  NUMBER  TWO  IN  DATA  FILE  UNITY  INPUT 

uydata ( : , 2) =ones (ndp, 1) ; 

% COLUMN  NUMBER  THREE  IN  DATA  FILE  ANGLE  OF  ATTACK 

col3= [ ' uydata ( : , 3) =simdata ( : , 2) ; ' ] ; eval (col3) ; 

% COLUMN  NUMBER  FOUR  IN  DATA  FILE  PITCH  RATE 

col4= [ ' uydata ( : , 4) =simdata ( : , 3) ; ' ] ; eval (col4) ; 

% COLUMN  NUMBER  FIVE  IN  DATA  FILE  THETA 

col5= [ 'uydata ( : , 5) =simdata ( : , 4) ; ' ] ; eval (col5) ; 

% COLUMN  NUMBER  SIX  IN  DATA  FILE  NORMAL  ACC 

col6= [ 'uydata ( : , 6) =simdata ( : , 5) ; ' ] ;eval (col6) ; 

% SENSOR  PLACEMENT  CORRECTIONS  FOR  AOA  AND  NORMAL  ACC 

uydata ( : , 3) =uydata ( : , 3) + (Xap*uydata ( : , 4) /vtrue) ; 
for  i=l : ndp; 

qsqr (i)= [uydata (i, 4) ] A2; 

q2=qsqr' ; 
end 

uydata ( : , 6) =uydata ( : , 6) - (Zan*q2/32 . 17) ; 
qdot=[diff (uydata (: ,4) ) * (1/dt) ;0] ; 
uydata ( : , 6) =uydata ( : , 6) - (Xan*qdot/32 .  17)  ; 
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,1)  =  DELTA  E  (INPUT) 

,3)  =  AOA 

,4)  =  PITCH  RATE  (Q) 

,5)  =  THETA 

, 6 )  =  NORMAL  ACC 


% INITIAL  CONDITIONS 

del=uydata (1,1) ;all=uydata (1,3) ;ql=uydata (1,4) ; 

thl=uydata (1,5) ;anl=uydata (1,6); 

% ADDITIONAL  INPUTS  TO  MMLE  FOLLOW 

p2snam='npsp2ss3' ;  %  MACRO  NAME  FOR  P2SS  FUNCTION 

pO=pref;  %-  INITIAL  PARAMETER  ESTIMATES  INPUT  DURING  NPSINIT3.M 

% —  CHECK  IF  THE  WEIGHTING  FUNCTION  IS  TO  BE  USED  FOR  INITIAL 

VALUES 

disp ('DO  YOU  WANT  TO  WEIGHT  THE  INITIAL  ESTIMATES?  INPUT  1=Y  0=NO 

') 

input ( '  ' ) ; 

if  ans==l 

dispC  Input  Weighting  Row  Vector  length  1x5') 

disp ('Use  brackets-  ex.  [.111. 11]   &  lower  #  higher  weight') 

rmsO=input ('  ' ) ; 

end 

pidq=[l];% IDENTIFY  WHICH  PARAMETERS  ARE  TO  BE  IDENTIFIED 

pidm=[l:5];% IN  THE  QUADRATIC,  MARQUARDT,  AND  FINAL  STAGES. 

pidf=[l:5] ;%-  pidq,  pidm,  pidf  MUST  BE  VALID  EVEN  IF  NOT  USED 

opt=[0  5  5  10  .02  .10  .001  1]; 

%-  DEFAULT  ITERATIONS  AND  CONVERGENCE  CRITERIA,  IF  NOISE  FREE 

OPT4=0 

gg0=eye (4) * ( . 001) ; % INNOVATIONS  COVARIANCE  MATRIX 

pert =le-4;% -PERTURBATION  USED  FOR  NUMERICAL  GRADIENT  CALCULATION 

linesearch=l;% USE  LINESEARCH  TO  HELP  PROC .  NOISE  CONVERGENCE 

!cls 

iranle    % CALL   MAIN   MMLE   MACRO   FROM   TOOL   BOX 

% PERFORM  FINAL  CALCULATIONS 

cla=pfin(l) ; cma=pf in (2) ; cmq=pf in (3) ; clde=pf in (4) ; cmde=pf in (5) ; 

dispC  ') 

dispC  '  ) 

deriv=[cla  cma  cmq  clde  cmde] ; 

dispC  MMLE    STABILITY  &  CONTROL  DERIVATIVES       ') 

dispC      CLA        CMA        CMQ        CLDE       CMDE') 

disp (deriv) 

dispC  INITIAL  INPUT  DERIVATIVES      ') 

disp (pref ) 

pause 

mleplot3 

diary  off 

%! print  npsmmle3.log; 

% END  NPSMMLE3.M 
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2.   NPSINIT3.M 

clear; 

%  MACRO  FILE  NAME  >=======  NPSINIT3.M  =======< 

%  Date:  31  Dec  91 

%  INITIAL  SETUP  MACRO  FOR  RUNNING  THE  PARAMETER 

%  ESTIMATION  TOOLBOX  IN  MATLAB  FOR  SIMPLIFIED  LONGITUDINAL 

%  SHORT  PERIOD  STABILITY  AND  CONTROL  DERIVATIVES 

%  THIS  MACRO  'GETS'  CONSTANTS  AND  DATA  FILE 

% LOAD  DATA  FILE 

dispC  DATA  FILE  MUST  CONTAIN  DATA  IN  COLUMN  MATRIX   ') 

dispC  MATRIX  NAME  SIMDATA:  N  DATA  PTS  X  5  COLUMNS  ') 

dispC  ELEVATOR/ALPHA/THETA/PITCH  RATE/NORMAL  ACC . ' ) 

dispC  ') 

dispC  DATA  FILE  NAME — MUST  EXIST  WITH  A   .mat   EXTENSION.  '); 

data=input ('ENTER  DATA  FILE  NAME  (  ==>  WITH  OUT  <==  .MAT 

EXTENSION) ?   ' , ' s' ) ; 

if  exist ('dt' ) ==0, dt=input ('DELTA  T  BETWEEN  DATA  POINTS?  ' ) ; end 

ldc=['load  ' , data, ' .mat; '] ; 

eval(ldc);% — EXECUTES  LOAD  COMMAND 

ndp=length (simdata) ; 

t=[0:ndp-l] *dt; 

% INPUT  REQUIRED  CONSTANTS IF  NOT  IN  DATA  FILE 

if  exist ('sref ) ==0, sref=input ('REFERENCE  AREA  (S)  IN  SQUARE  FEET? 

' ) ; end 

if  exist ('cbar' ) ==0, cbar=input ('MEAN  AERODYNAMIC  CHORD  IN  FEET? 

' ) ; end 

if  exist ('gw' ) ==0, gw=input ('AIRCRAFT  GROSS  WEIGHT  IN  POUNDS?  ');end 

if  exist (' iyy' )==0, iyy=input ('MOMENT  OF  INERTIA  (IYY)  IN  SLUG-FTA2? 

' ) ; end 

if  exist ('Xap' )==0,Xap=input ('X-DIST  FROM  eg  TO  AOA  PROBE  (FT 

+FWD) ' ) ;end 

if  exist ('Zan' ) ==0, Zan=input ('Z-DIST  FROM  eg  TO  NORMAL  ACCEL  (FT 

+DWN) ' ) ;end 

if  exist ('Xan' ) ==0, Xan=input ('X-DIST  FROM  eg  TO  NORMAL  ACCEL  (FT 

+FWD ) ' ) ; end 

%sdc=['save   ',data,'   simdata  dt  sref  cbar  gw  iyy  Xap  Zan 

Xan' ] ;eval (sdc) ; 

vtrue=input ('AIRSPEED  IN  FEET  PER  SECOND?  '); 

altft=input ('AIRCRAFT  ALTITUDE  IN  FEET?   '); 

oat=input ('OUTSIDE  AIR  TEMPERATURE?     '); 

% CALCULATE  CONSTANTS  DENSITY  AND  DYNAMIC  PRESSURE 

rho=.0023769*exp( (-1*32 . 174*altft) / (1716*  (oat+460) ) ) ; 
qbar= . 5*rho*vtrue*vtrue; 
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% INPUT  INITIAL  ESTIMATES  FOR  STABILITY 

% AND  CONTROL  DERIVATIVES 

pref (l)=input ('CL_ALPHA  ESTIMATE  (1/RAD) ?  '); 

pref (2) =input ('CM_ALPHA  ESTIMATE  (1/RAD)?  '); 

pref  (3)  =input  ('CM_Q_  ESTIMATE  (1/RAD)?  '); 

pref (4)=input ( ' CL_DE  ESTIMATE  (1/RAD)?  '); 

pref (5)=input ( ' CM_DE  ESTIMATE  (1/RAD)?  '); 

.'erase  npsinit3  .mat; 

save  npsinit3 

% END  NPSINIT3.M 
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3.   NPSP2SS3.M 


function 
%   MACRO 

% 

% 


[a, phi, gam, c, d, q, xO, dt, rowinq,b] =npsp2ss3 (p) 
FILE  NAME     >=======  NPSP2SS3.M  =======< 

Date:  31  Jan  92 


MACRO  TO  EST.  FUNCTION  FOR  TRANSFORMING 
MODEL  PARAMETERS  INTO  STATE  SPACE  EQUATIONS 


P2SS  FUNCTION  FOR  NPSMMLE3.M 


P(D 
P(2) 
P(3) 
P(4) 
P(5) 


CL 
CM~ 

cm" 

CL~ 
CM 


ALPHA 
"ALPHA 

Q 

DE 
DE 


STABILITY  AND  CONTROL 
PARAMETERS 


PERFORM  INITIAL 

constl= (-l*qbar*sref ) / (cos (all) *gw*vtrue/32 . 17) ; 

const 2=qbar*sref*cbar/iyy; 

const3=const2*cbar/ (2*vtrue) ; 

const4=32 . 17*cos (thl) / (vtrue*cos (all)  )  / 

const5= (qbar*sref ) /gw; 

const6=-l* (constl*p (1) *all+ql+constl*p (4) *del+const4) ; 

const7=-l*  (const2*p (2) *all+const3*p (3) *ql+const2*p (5) *del) ; 

const8= (-1* (const5*p (1) *all+const5*p (4) *del) ) +anl; 

% STABILITY  DERIVATIVES 

a= [constl*p (1)  1 

const2*p(2)  const3*p(3) 

0  1 

% 

b= [constl*p (4)      (const4+const6) ; 
const2*p(5)  const7; 

0  <-l*ql)]; 


CALCULATIONS 


0; 

0; 
0]; 
CONTROL  DERIVATIVES 


c=[l 
0 
0 

const5*p (1) 

% 

d=[0 
0 
0 
const5*p (4) 


-MEASUREMENT 

0, 

0, 

1, 
0] 
FEED 


MATRIX 


THROUGH  MATRIX 


0 

0 

0 

const8] 


% STATE  NOISE  COVARIANCE 

q=eye (a) *le-4;% —  Q  IS  THE  SAME  SIZE  AS  a 

% WITH  Q*Q'  POS.  DEFINITE 

% ROWS  IN  Q  IN  WHICH  PARAMETERS  OCCUR,  A  VECTOR 

%  SAME  DIMENSION  AS  p 
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rowinq= [0  0  0  0  0 ] ; 

% INITIAL  STATE  VECTOR 

x0=[all  ql  thl] ; 

% DISCRETIZE 

%  *****NEED  TO  EDIT  dt  (below)  FOR  THE  DELTA  T  OF  THE  DATA  ***** 
dt= . 1 / 

[phi, gam] =c2d (a,b, dt) / 
% END  NPSP2SS3.M 
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MLEPL0T3.M 


%  MACRO  FILE  NAME         >=======  MLEPLOT3.M  =======< 

%  Date:  31  Jan  92 

%  MACRO  TO  PLOT  DATA  FROM  NPSMMLE3 

constl= (-l*qbar*sref ) / (cos (all) *gw*vtrue/32 . 17) ; 

const2=qbar*sref *cbar/iyy; 

const3=const2*cbar/ (2*vtrue) ; 

const4=32 . 17*cos (thl) / (vtrue*cos (all)  )  ; 

const5= (qbar*sref ) /gw; 

const 6=-l* (constl*pf in (1) *all+ql+constl*pf in (4) *del+const4) ; 

const7=-l* (const2*pf in (2) *all+const3*pf in (3) *ql+const2*pf in (5) *del) 


const 8= (-1* (const5*pf in (1) *all+const5*pf in (4) *del) ) +anl; 

% STABILITY  DERIVATIVES 

a=[constl*pfin(l)  1  0 

const2*pfin(2)  const3*pf in (3)  0; 

0  10] 

% CONTROL  DERIVATIVES 

b= [constl*pf in (4)      (const4+const6) ; 

const2*pf in  (5)  const7; 

0  (-l*ql)]; 

% MEASUREMENT  MATRIX 

c=[l  0  0, 

0  10 

0  0  1 

const5*pfin  (1)              0  0] 

% FEED  THROUGH  MATRIX 


const8] 


d=[0 

0 

0 

const5*pf in (4) 
rtdc=(180/pi) ; 

% OUT2  =  OUTPUT  VECTOR OUT3  =  STATE  VECTOR 

[OUT2,OUT3]=lsim(a,b,c,d,uydata(: , 1:2) ,t,x0) ; 

% PLOTS  TO  MONITOR  AND  STORED  IN  META  FILE 

if  exist ('typac' ) ==0, typac=input (' INPUT  THE  AIRCRAFT  TYPE  ? 


'  ,'3 


) 


! erase 
! erase 
! erase 
! erase 
! erase 


a 
a 
a 
a 


end 

a : \plots\outputg.met; 

\plot s\outputth . met ; 

\plots\outputde . met ; 

\plot s\outaoa . met ; 

\plots\outputq.met ; 

% ELEVATOR  VS  TIME 

hold  off ;plot (t,rtdc*uydata( : , 1) , ' -t' ) ; 

xlabel('Time  (seconds) ') /ylabel (' Elevator  Input  (degrees)'); 

ans=[' title  ('  '  '  ,  typac,  '  ELEVATOR  INPUT  VS  TIME  ");']; 

eval (ans) 

pause 

meta  A:\plots\OUTPUTde 
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% AOA  (OBSERVED  AND  ESTIMATED)  VS  TIME 

plot (t, rtdc*uydata ( : , 3) , ' *r' ) ;hold  on; 

xlabel (' Time  (seconds) ') ;y label ( 'AOA  (degrees)'); 

ans=[' title  (' ' ' , typac, '   ESTIMATED  AND  MEASURED  AOA  RESPONSE'');']; 

eval (ans) 

text  (  . 6,  . 85, '  *   Measured  Data  Points  ' , ' sc' ) ;pause; 

plot (t,rtdc*OUT2 (: , 1) , ' og' ) ; 

text  (  . 6,  . 80, ' o   Estimated  Response  ' , ' sc' ) ;pause; 

pause 

meta  A:\plots\outAOA 

% Q  (OBSERVED  AND  ESTIMATED)  VS  TIME 

hold  off;plot (t, rtdc*uydata(:,4) , ' *r' ) ;hold  on; 

xlabel ('Time  (seconds) '); ylabel (' Pitch  Rate,  Q,  (deg/sec) ' ) ; 

ans= [ 'title (''' ,typac, '   ESTIMATED  AND  MEASURED  q  RESPONSE'');']; 

eval (ans) 

text  ( . 6,  . 85, ' *   Measured  Data  Points  ' , ' sc' ) ;pause; 

plot (t, rtdc*0UT2 ( : , 2) , ' og' ) ; 

text  ( . 6,  . 80, ' o   Estimated  Response  ',' sc' ); pause; 

pause 

meta  A:\plots\OUTPUTQ 

% THEATA  (OBSERVED  AND  ESTIMATED)  VS  TIME 

hold  of f ;plot (t, rtdc*uydata ( : , 5) , ' *r' ) ;hold  on; 

xlabel ('Time  (seconds) '); ylabel (' Pitch  Angle,  Theta,  (deg) ' ) ; 

ans=[' title (''' ,typac, '   ESTIMATED  AND  MEASURED  THETA 

RESPONSE' ');']; 

eval (ans) 

text  (  . 6,  . 85, ' *   Measured  Data  Points  ' , ' sc' ) ;pause; 

plot (t,rtdc*OUT2 ( : ,3) ,  '  og'  )  ; 

text  (  . 6,  . 80, ' o   Estimated  Response  ' , ' sc' ) /pause; 

pause 

meta  A:\plots\OUTPUTTH 

% ACCELERATION  (OBSERVED  AND  ESTIMATED)  VS  TIME 

hold  of f ;plot (t, uydata ( : , 6) , ' *r' ) ;hold  on; 

xlabel ('Time  (seconds) ') ;ylabel ( 'Acceleration,  G  '); 

ans= [' title (''' ,typac, '   ESTIMATED  AND  MEASURED  G  RESPONSE'');']; 

eval (ans) 

text  (  . 6,  . 85,  '  *   Measured  Data  Points  ' , ' sc' ) ;pause; 

plot(t,OUT2(:,4)  ,  '  og'  )  ; 

text ( . 6,  . 80,  '  o   Estimated  Response  ' , ' sc' ) ;pause; 

pause 

meta  A:\plots\OUTPUTG 

hold  off; 

% END  MLEPLOT3  .M 
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F.      ACTUAL  LATERAL-DIRECTIONAL  MMLE 

1.   NPSMMLE4  M 


MACRO  NAME 


% 

% 

clear; 

!  erase  npsnunle4.log; 

diary  npsmmle4.log 


>======  NPSMMLE4.M  =======< 

Date:  5   Feb   92 


STABILITY  PARAMETER  IDENTIFICATION  MACRO 
FOR  LATERAL-DIRECTIONAL 
STABILITY  AND  CONTROL  DERIVATIVES 


FILE' ) 
') 
') 


disp (ans) 
disp('NPS 
disp  (' 
disp  ( ' 
dispC  ') 
disp (ans) 

npsinit4  % RUN  INITIALIZATION  MACRO 

format  compact , clc; 

load  npsinit4; 

global  sref  bbar  gw  ixx  ixz  izz  vtrue  qbar  dt  betal  rolll  yawl 

ranglel  ayl  dal  drl; 

% INPUT  FLIGHT  TEST  DATA 

uydata=zeros (ndp, 8) ;%  ESTABLISH  DATA  MATRIX  FOR  MMLE.M 


%  UYDATA( 

(P) 

%  UYDATA( 

%  UYDATA( 

(phi) 

%  UYDATA( 

(ay) 

% 

coll  = 

% 

col2  = 


1)  =  DELTA  Aileron  (INPUT) 


UYDATA ( : , 5 )  =  ROLL  RATE 


,2) 
,3) 


=  DELTA 
=  UNITY 


Rudder 
INPUT 


(INPUT) 


UYDATA ( 
UYDATA ( 


,6) 
,7) 


=  YAW  RATE  (r) 
=  ROLL  ANGLE 


4)  =  BETA  (SIDE  SLIP; 


UYDATA ( : , 8 )  =  LATERAL  G 


NUMBER  ONE  IN 
: , 1) =simdata ( 


COLUMN 
' uydata ( 

COLUMN  NUMBER  TWO  IN 
' uydata (  :  ,  2 ) =simdata ( 
uydata ( : , 3) =ones (ndp, 1) ; % 

% 

col4  = 


DATA  FILE  AILERON  INPUT 
,1) ;'];eval(coll) ; 
DATA  FILE  RUDDER  INPUT 
,2) ; ' ] ;eval (col2) ; 
UNITY  INPUT 
COLUMN  NUMBER  FOUR  IN 
' uydata ( : , 4) =simdata ( 


DATA  FILE  BETA  (beta) 
,3);'];eval(col4) ; 


% COLUMN  NUMBER  FIVE  IN  DATA  FILE  ROLL  RATE  (p) 

col5= [ ' uydata ( : , 5) =simdata ( : , 4) ; ' ] ; eval (col5) / 

% COLUMN  NUMBER  SIX  IN  DATA  FILE  YAW  RATE  (r) 

col 6= [ ' uydata ( : , 6) =simdata ( : ,5) ; ' ] ; eval (col 6) / 

% COLUMN  NUMBER  SEVEN  IN  DATA  FILE  ROLL  ANGLE  (phi) 

col 7= [ 'uydata ( : , 7) =simdata ( : ,  6) ; ' ] ; eval (col7) ; 

% COLUMN  NUMBER  EIGHT  IN  DATA  FILE  LATERAL  G  (ay) 

col8= [ ' uydata ( : , 8) =simdata ( : , 7) ;  '  ] ;eval (col8)  ; 

% SENSOR  PLACEMENT  CORRECTIONS  FOR  beta  AND  ay 

uydata ( : ,  4) =uydata ( : ,  4) - (Xbp/vtrue) *uydata ( : ,  6) ; 

rdot=[diff (uydata (: , 6) ) * (1/dt) ;0] ; 

pdot=[diff (uydata <:, 5) )* (1/dt) ;0] ; 
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uydata ( : ,  8) =uydata ( : , 8) - (rdot*Xay/32 . 17) + (pdot*Zay/32 . 17) ; 

% INITIAL  CONDITIONS  FOR  da,  dr,  BETA,  p,  r,  phi,  ay 

betal=uydata (1,4) ;rolll=uydata (1,5) ;yawl=uydata (1,6); 
ranglel=uydata (1,7) ; ayl=uydata (1,8) ; 
dal=uydata (1,1) ;drl=uydata (1,2) ; 

% ADDITIONAL  INPUTS  TO  MMLE  FOLLOW 

p2snam='npsp2ss4' ;% MACRO  NAME  FOR  P2SS  FUNCTION 

pO=pref/% INITIAL  PARAMETER  ESTIMATES 

% — CHECK  IF  WEIGHTING  FUNCTION  DESIRED 

disp('Do  you  want  to  WEIGHT  the  initial  estimates?  INPUT  1=YES 

0=N0' ) ; 

ans=input ('  '  )  ; 
if  ans==l, 

disp ('Input  WEIGHTING  ROW  MATRIX:   1  X  12'') 

disp ('USE  BRACKETS-  ex.  [.1  .1  1  1  1  10  1  .1  1  1  .1  1]'); 

disp ('NOTE  ***  The  LOWER  the  #  the  HIGHER  the  WEIGHTING!'); 

rmsO=input ('  '  )  ; 
end; clc; 

pidq=[l];% IDENTIFY  WHICH  PARAMETERS  ARE  TO  BE  IDENTIFIED 

pidm=[l:12]  ;% IN  THE  QUADRATIC,  MARQUARDT,  AND  FINAL 

STAGES . 

pidf=[l:12] ;% pidq,  pidm,  pidf  MUST  BE  VALID  EVEN  IF  NOT  USED 

opt=[0  5  5  10  .02  .10  .001  1];%-  DEFAULT  ITERATIONS  AND  CONVERGENCE 
%  CRITERIA,  IF  NOISE  FREE  OPT4=0 

gg0=eye  (5)  *  (  .  01)  ;% INNOVATIONS  COVARIANCE  MATRIX 

pert=le-4;% —  PERTURBATION  USED  FOR  NUMERICAL  GRADIENT  CALCULATION 

linesearch=l;% USE  LINESEARCH  TO  HELP  PROC .  NOISE  CONVERGENCE 

mmle% CALL  MAIN  MMLE  MACRO  FROM  TOOL  BOX 

% PERFORM  FINAL  CALCULATIONS 

CYb=pfin(l) ;Clb=pfin(2) ;CNb=pfin(3) ;Clp=pfin(4) ; 
CNp=pfin(5) ;Clr=pfin(6) ;CNr=pfin(7) ;Clda=pf in (8) ; 
CNda=pfin(9) ;CYdr=pf in (10) ;Cldr=pf in (11) ;CNdr=pf in (12) ; 
deriv=[CYb        Clb         CNb        Clp        CNp        Clr; 

CNr        Clda        CNda       CYdr       Cldr       CNdr] ; 
disp('  MMLE    STABILITY  &  CONTROL  DERIVATIVES        ') 

disp('CY_b        Cl_b        CN_b        Cl_p        CN_p      Cl_r  ') 
disp (deriv (1, : ) ) 

disp('CN_r       Clda       CNda       CYdr       Cldr      CNdr') 
disp (deriv (2, : ) ) 

dispC  INITIAL  INPUT  DERIVATIVES      ') 

disp (pref (1:6)); 
disp (pref (7:12) ) ; 
pause;mleplot4 
diary  off 

%! print  npsmmle4.log; 
% END  NPSMMLE4.M 
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NPSINIT4.M 


clear; 

%  MACRO  FILE  NAME  >=======  NPSINIT4.M  =======< 

%  Date:  5   Feb  92 

%  INITIAL  SETUP  MACRO  FOR  RUNNING  THE  PARAMETER 

%  ESTIMATION  TOOLBOX  IN  MATLAB  FOR  LATERAL-DIRECTIONAL 

%  STABILITY  AND  CONTROL  DERIVATIVES 

%  THIS  MACRO  'GETS'  CONSTANTS  AND  DATA  FILE 

disp('  ') 

dispC   DATA  FILE  NAME — MUST  EXIST  WITH  A   .mat   EXTENSION.  '); 

dispC     DATA  FILE  MUST  CONTAIN  DATA  IN  COLUMN  MATRIX'); 

dispC    MATRIX  NAME  SIMDATA:   N  (data  pts)  X  7  COLUMNS'); 

disp ('AILERON /RUDDER/ BETA/ ROLL  RATE/YAW  RATE/ROLL  ANGLE/LAT  G' ) ; 

data=input ('ENTER  DATA  FILE  NAME  (  ==>  WITH  OUT  <==  .MAT 

EXTENSION)?   ','s'); 

ldc= [ ' load  ' , data,  '  .mat; ' ] ; eval (ldc) ; 

ndp=length (simdata) ; 


if  exist 


t=[0:ndp-l]*dt; 

% 


if  exist 
' ) ; end 
if  exist 
if  exist 
if  exist 
' ) ; end 
if  exist 
' ) ; end 
if  exist 
' ) ; end 
if  exist 


+FWD) ' ) ;end 


if  exist 


' dt ' ) ==0 , dt=input ( ' INPUT  dt  BETWEEN  DATA  POINTS .   ' ) ; end 


INPUT  REQUIRED  CONSTANTS 

'sref )==0,sref=input ('REFERENCE  AREA  (S)  IN  SQUARE  FEET? 

'bbar' ) ==0,bbar=input ('WINGSPAN  IN  FEET?  ');end 

'gw' )==0,gw=input ('AIRCRAFT  GROSS  WEIGHT  IN  POUNDS?  ' ) ; end 

' ixx' )==0, ixx=input ('MOMENT  OF  INERTIA  (Ixx)  IN  SLUG-FT^2? 

' ixz' ) ==0, ixz=input ('MOMENT  OF  INERTIA  (Ixz)  IN  SLUG-FTA2? 

' izz' ) ==0, izz=input ('MOMENT  OF  INERTIA  (Izz)  IN  SLUG-FTA2? 


'Xbp' )==0,Xbp=input ('X-DIST  FROM  eg  TO  BETA  PROBE  (FT 


' Zay' )==0, Zay=input ('Z-DIST  FROM  eg  TO  LATERAL  ACC .  (FT 


+DWN) ' ) ;end 

if  exist ('Xay' ) ==0,Xay=input ('X-DIST  FROM  eg  TO  LATERAL  ACC.  (FT 

+FWD ) ' ) ; end 

%sdc=['save   ',data,'   simdata  dt  sref  bbar  gw  ixx  ixz  izz  Xbp  Zay 

Xay' ]  ; 

%eval (sdc) ; 

vtrue=input ('AIRSPEED  IN  FEET  PER  SECOND?  '); 

altft=input ('AIRCRAFT  ALTITUDE  IN  FEET?   '); 

oat=input ('OUTSIDE  AIR  TEMPERATURE?     '); 

rho=.0023769*exp( (-1*32 . 174*altft) / (1716* (oat+460) ) ) ; 

qbar= . 5*rho*vtrue*vtrue; 

% INPUT  INITIAL  ESTIMATES  FOR  STABILITY 

% AND  CONTROL  DERIVATIVES  TO  START  MMLE  PROGRAM 

pref  (l)=-.6;%  input ( ' CY_beta  FROM  WIND  TUNNEL  (1/RAD) ?  '); 
pref (2)=-.15;%  input ( ' Cl_beta  FROM  WIND  TUNNEL  (1/RAD)?  '); 
pref (3)=. 20;%  input (' CN  beta  FROM  WIND  TUNNEL  (1/RAD)?  '); 
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pref (4)=-.35;%  input (' Cl_p  FROM  WIND  TUNNEL  (1/RAD) ?  ') 

pref (5)=-.05;%  input (' CH_p  FROM  WIND  TUNNEL  (1/RAD)?  ') 

pref (6)=. 15;%  input('Cl_r  FROM  WIND  TUNNEL  (1/RAD)?  ') 

pref (7)=-. 2;%  input (' CN_r  FROM  WIND  TUNNEL  (1/RAD)?  ') 

pref (8)=. 05;%  input ( ' Cl_da  FROM  WIND  TUNNEL  (1/RAD)?  ') 

pref (9)=-. 001;%  input (' CN_da  FROM  WIND  TUNNEL  (1/RAD)?  '); 

pref (10)=.175;%input ('CY_dr  FROM  WIND  TUNNEL  (1/RAD)?  '); 

pref (ll)=.02;%input ( ' Cl_dr  FROM  WIND  TUNNEL  (1/RAD)?  '); 

pref (12)=-.075;%input ( ' CN_dr  FROM  WIND  TUNNEL  (1/RAD)?  '); 
! erase  npsinit4 .mat; 
save  npsinit4 

% END  NPSINIT4.M 
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3.   NPSP2SS4.M 


function  [a, phi, gam, c, d, q, xO, dt, rowinq, b] =npsp2ss4 (p) 
%   MACRO  FILE  NAME     >=======  NPSP2SS4.M  =======< 

%  Date:  5   Feb  92 


% 
% 

% 

%  P2SS  FUNCTION  FOR  NPSMMLE4.M 


MACRO  TO  ESTABLISH  FUNCTION  FOR  TRANSFORMING 
MODEL  PARAMETERS  INTO  STATE  SPACE  EQUATIONS 


(1) 
(2) 
(3) 
(4) 
(5) 
(6) 


CY_beta 

Cl_beta 

CN_beta 

Cl_p 

CNjp 

CI  r 


STABILITY  AND  CONTROL 


PARAMETERS 


p(7)  =  CN_r 
p(8)  -  Cl_da 
p(9)  =  CN_da 
p(10)  =  CY_dr 
p(ll)  =  Cl_dr 
p(12)  =  Cl_dr 
CALCULATIONS 


% PERFORM  INITIAL 

constl= (qbar*sref ) / (gw/32 . 17) ; 

const2=qbar*sref *bbar ; const2a=const2*bbar/ (2*vtrue) 

const3= (qbar*sref ) /gw; 

% INERTIAL  MATRIX 


In= 


[1 
0 
0 
0 


0 

ixx 
-ixz 
0 


0 
—ixz 
izz 
0 


0; 
0; 
0; 
1] 


% PLANT 

an= [constl*p (1) /vtrue  0 

const2*p(2)  const2a*p(4) 
const2*p(3)  const2a*p(5) 
0  1 

% 

a=inv (In) *an; 

% 

bn=[0       constl*p (10) /vtrue 
const2*p(8)    const2*p(ll) 


const2*p(9) 

0 
% 
b=inv(In) *bn; 


const2*p(12) 
0 


-1 

(32 

17/vt. 

cue) 

const2a 

*P(6) 

0 

const2a 

*p(7) 

0 

0 

0] 

0 
0 
0, 
0]/ 


c=[l 

0 

0 

0 

const3*p  (1) 
% 
d=[0 

0 

0 


0 

0 

o, 

1 

0 

o, 

0 

1 

o( 

0 

0 

1, 

0 

0 

0], 

0 

betal; 

0 

rolll; 

0 

] 

I  awl ; 

98 


0  0  ranglel; 

0  const3*p(10)       ayl]; 

% 

% STATE  NOISE  COVARIANCE 

q=eye(a)  *le-4;% Q  IS  THE  SAME  SIZE  AS  a 

%  WITH  Q*Q'  POS.  DEFINITE! 

% ROWS  IN  Q  IN  WHICH  PARAMETERS  OCCUR,  A  VECTOR 

%  SAME  DIMENSION  AS  p 

rowinq=0*p; 

% INITIAL  STATE  VECTOR 

xO=[betal  rolll  yawl  ranglel]; 

% DISCRETIZE 

%*****N0TE*****  CHANGE  dt  TO  THE  ACTUAL  DATA  VALUE  **** 

dt=.05; 

[phi, gam] =c2d (a,b, dt) ; 

% END  NPSP2SS4.M 
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MLEPL0T4 . M 


%  MACRO  FILE  NAME 

% 

% 

% 

% 


>=======  MLEPLOT4.M  =======< 

Date:  5   Feb   92 

MACRO  TO  PLOT  DATA  FROM  NPSMMLE4 
GENERATE  THE  PREDICTED  DATA 
CALCULATE  STABILITY  AND  CONTROL  MATRIX 


mass=gw/32 . 17;rtdc= (180 /pi) ; const 1= (qbar*sref ) /mass; 
const2=qbar*sref *bbar;const2a=const2*bbar/ (2*vtrue) ; 
const 3=qbar*sref/gw; 

% 


In=[l    0    0    0; 

0   ixx  —ixz  0; 

0  -ixz   izz  0; 

0    0    0    1  ]  ; 
an= [constl*pf in (1) /vt 

const2*pfin(2) 

const2*pf in (3) 

0 
a=inv (In) *an; 

2- 

rue        0 

const2a*pf in (4)   cons 
const2a*pf in (5)   cons 
1 

-1 
t2a 
t2a 
0 

*pf. 
*pf. 

32. 
Ln(6) 
Ln(7) 

17 /vtrue; 
0; 
0; 
0]/ 

bn=[0 

const2*pfin (8) 
const2*pfin (9) 
0 

b=inv (In) *bn; 

% 

c=[l 
0 
0 
0 
const3*pf in (1) 

% 

d=[0          0 
0           0 
0           0 
0           0 
0           consi 

constl*pf in (10) /vtrue 

const2*pfin(ll) 

const2*pfin (12) 

0 

0, 
0, 

o, 

0], 

0 
1 
0 
0 
0 

0 
0 
1 
0 
0 

0, 
0, 
0, 
1, 
0], 

t3* 

pfin(10) 

betal; 
rolll; 

yawl  ; 
ranglel; 

ayl]  ; 

%  OUT2  =  OUTPUTS  — 



OUT3  = 

STATE  1 

7ECTOR 

[OUT2,OUT3]=lsim(a,b,c,d,uydata ( : , 1 


3) ,t,x0) ; 
AND  STORE 


TO  META  FILE 


%     PLOTS  FOR  VIEWING  ON  MONITOR 
.'erase  a  :  \plots\*  .met 

subplot (211) ;plot (t,rtdc*uydata ( : , 1) ) ; 
xlabel('Time  (seconds) ') ;ylabel ( 'Aileron  Input 
ans=['  title  ('"  ,typac,  '  AILERON  INPUT  VS  TIME 
subplot (212) ;plot (t , rtdc*uydata ( : ,2) ) ; 
xlabel('Time  (seconds) ') ;ylabel (' Rudder  Input  (degrees) 


(degrees) ' ) ; 
' ) ; ' ] ;eval (ans) 


); 
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ans=[' title  (' ' ' , typac, '  RUDDER  INPUT  VS  TIME  "  ) ; ' ] ;eval (ans) 

pause;  %meta  A:\plots\lNPUTDAR 

%   Beta  vs  Time 

subplot (111) /plot (t, rtdc*uydata ( : , 4) , ' *r' ) ;hold  on; 

xlabel('Time  (seconds) ' ) ;ylabel ('Beta  (degrees)'); 

ans=  [' title  ('",  typac,  '   ESTIMATED  AND  MEASURED  Beta 

RESPONSE' ');']; 

eval (ans) 

text ( . 6, . 85, ' *   Measured  Data  Points  'f 'sc') /pause; 

plot (t, rtdc*0UT2 ( : , 1) , ' og' ) ; 

text  ( . 6,  . 80, ' o   Estimated  Response  ' , ' sc' ) ;pause; 

pause;  %meta  A:\plots\outbeta 

%  Roll  rate  vs  Time 

hold  of f ;plot (t, rtdc*uydata ( : , 5) , ' *r' ) ;hold  on; 

xlabel('Time  (seconds) ') ;ylabel (' Roll  Rate,  p,  (deg/sec)')/ 

ans= [' title (''' , typac, '   ESTIMATED  AND  MEASURED  p  RESPONSE'');' 

eval (ans) 

text  (  . 6,  . 85, ' *   Measured  Data  Points  ' , ' sc' ) ;pause; 

plot (t, rtdc*OUT2 ( : , 2) , ' og' ) / 

text ( . 6, . 8  0, ' o   Estimated  Response  ' , ' sc' ) ; pause; 

pause; %meta  A:\plots\OUTP 

%   Yaw  Rate  vs  Time 

hold  of f ;plot (t, rtdc*uydata ( : , 6) , ' *r' ) ;hold  on; 

xlabel('Time  (seconds) ') ;ylabel (' Yaw  Rate,  r,  (deg/sec)')/ 

ans= [ 'title ("' , typac, '   ESTIMATED  AND  MEASURED  r  RESPONSE'');' 

eval (ans) 

text ( . 6, . 85, ' *   Measured  Data  Points  ' , ' sc' ) ;pause; 

plot (t,rtdc*OUT2 ( : , 3) , ' og' ) ; 

text ( . 6, . 80, ' o   Estimated  Response  ' , ' sc' ) ;pause; 

pause; %meta  A:\plots\OUTr 

%  Bank  Angle  vs  Time 

hold  off /plot (t,rtdc*uydata(: ,7) , '  *r' ) ;hold  on; 

xlabel('Time  (seconds) ') /ylabel (' Bank  Angle,  phi,  (deg) ' ) / 

ans=[' title (' ' ' , typac, '   ESTIMATED  AND  MEASURED  phi  RESPONSE'') 

eval (ans) 

text ( . 6, . 85, ' *   Measured  Data  Points  ',' sc' ) /pause/ 

plot (t, rtdc*OUT2 ( : , 4) , ' og' ) / 

text ( . 6, . 80, ' o   Estimated  Response  ',' sc' ) /pause/ 

pause/  %meta  A:\plots\OUTphi 
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%  Lateral  G  vs  Time 

hold  off ;plot (t,uydata( : ,  8) , ' *r' ) ;hold  on; 

xlabel (' Time  (seconds) ') ;ylabel (' Lateral  G,  ay,  (G) ' ) ; 

ans=  ['  title  ('"  ,typac,  '   ESTIMATED  AND  MEASURED  Lateral 

RESPONSE' ');']; 

eval (ans) 

text ( . 6, . 85, ' *   Measured  Data  Points  ' , 'sc' ) /pause; 

plot (t , 0UT2 ( : , 5) , ' og' ) ; 

text  ( . 6,  . 80, ' o   Estimated  Response  ' , 'sc') ; pause; 

pause;  %meta  A:\plots\OUTlatg 

hold  off; 

% END  MLEPLOT4 .M 
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APPENDIX  B  OUTPUT 
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A-4D  SIMULATED  LONGITUDINAL  RESULTS 
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Figure   B-9   Navion   Pitch  Angle   Response 
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Figure  B-10  Navion  Normal  Acceleration  Response 
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NOTION    SIMULATED    LONGITUDINAL    RESULTS 
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Figure  B-14  UAV  Pitch  Angle  Response 
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Figure  B-15  UAV  Normal  Acceleration  Response 
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UAV  SIMULATED  LONGITUDINAL  RESULTS 

pid  p (pid)  pref  cramer  2fcramer  insens 

1.0000  5.1314  6.0000  0.0989     0.2913  0.0888 

2.0000  -0.9390  -0.6000  0.0244     0.0718  0.0177 

3.0000  -3.4725  -10.0000  1.2628     3.7180  0.6143 

4.0000  0.0473  0.5000  0.0248     0.0731  0.0232 

5.0000  -0.2188  -1.1000  0.0096     0.0283  0.0055 


MMLE    STABILITY  &  CONTROL  DERIVATIVES 

CLA         CMA         CMQ         CLDE  CMDE 

5.1314    -0.9390    -3.4725     0.0473  -0.2188 

INITIAL  INPUT  DERIVATIVES 

6.0000    -0.6000   -10.0000     0.5000  -1.1000 

"TRUTH  DERIVATIVES"  USED  TO  GENERATE  DATA 

5.0100    -0.6900   -14.7500     0.0760  -0.3333 
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Figure   B-19   A-4D   Bank  Angle   Response 
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Figure  B-20  A-4D  Roll  Rate  Response 
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Figure  B-21  A-4D  Lateral  Acceleration  Response 
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A-4D  SIMULATED  LATERAL -DIRECTIONAL  RESULTS 


pid 

p(pid) 

pref 

cramer 

2fcramer 

insens 

1.0000 

-0.9863 

-0.6000 

0.0102 

0.0273 

0.0085 

2.0000 

-0.1225 

-0.1500 

0.0011 

0.0030 

0.0007 

3.0000 

0.2512 

0.2000 

0.0011 

0.0030 

0.0008 

4.0000 

-0.2610 

-0.3500 

0.0030 

0.0079 

0.0013 

5.0000 

0.0143 

-0.0500 

0.0044 

0.0117 

0.0020 

6.0000 

0.1432 

0.1500 

0.0085 

0.0227 

0.0057 

7.0000 

-0.4134 

-0.2000 

0.0119 

0.0319 

0.0074 

8.0000 

0.0775 

0.0500 

0.0008 

0.0020 

0.0004 

9.0000 

0.0648 

-0.0010 

0.0012 

0.0032 

0.0006 

10.0000 

0.1818 

0.1750 

0.0072 

0.0194 

0.0068 

11.0000 

-0.1035 

0.0200 

0.0007 

0.0019 

0.0004 

12.0000 

0.0340 

-0.0750 

0.0010 

0.0026 

0.0006 

MMLB    STABILITY  6  CONTROL  DERIVATIVES 

CYJb        Cljb        CNJb        Cl_p        CN_p  Cl_r 

-0.9863    -0.1225     0.2512      -0.2610     0.0143  0.1432 

CN_r       Clda       CNda       CYdr       Cldr  CNdr 

-0.4134     0.0775     0.0648       0.1818    -0.1035  0.0340 

INITIAL  INPUT  DERIVATIVES 

-0.6000    -0.1500     0.2000      -0.3500    -0.0500  0.1500 

-0.2000     0.0500    -0.0010       0.1750     0.0200  -0.0750 

"TRUTH  DERIVATIVES"  USED  TO  GENERATE  DATA 

-0.9800    -0.1200     0.2500      -0.2600     0.0220  0.1400 

-0.3500     0.0800     0.0600       0.1700    -0.1050  0.0320 
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Figure  B-22  Navion  Aileron  and  Rudder  Inputs 
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Figure  B-23  Navion  Beta  Response 
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Figure  B-24  Navion  Yaw  Rate  Response 


NAVION 

ESTIMATED  AND  MEASURED  phi  RESPONSE 

3 
2 

*      Measured    Data    Points 
o     Estimated  R»«pon»« 

1 

—     "True  Response" 

•a 

0< 

-£. 

Q. 

-1 

f 

-2 

B 
o 
00 

-3 

-4 

-5 

-6 
< 

>                      2                     4 

6                    8                   10                  12                  14 

1 

6 

Time   (seconds) 

Figure  B-25  Navion  Bank  Angle  Response 
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NOTION  SIMULATED  LATERAL-DIRECTIONAL  RESULTS 

pid      p (pid)      pref     cramer    2fcramer  insens 

1.0000    -0.5630    -0.6000     0.0059     0.0171  0.0057 

2.0000    -0.0746    -0.1500     0.0017     0.0049  0.0006 

3.0000     0.0715     0.2000     0.0014     0.0040  0.0004 

4.0000    -0.4030    -0.3500     0.0086     0.0249  0.0019 

5.0000    -0.0602    -0.0500     0.0070     0.0201  0.0016 

6.0000     0.1029     0.1500     0.0055     0.0160  0.0028 

7.0000    -0.1189    -0.2000     0.0045     0.0129  0.0019 

8.0000     0.1307     0.0500     0.0028     0.0082  0.0009 

9.0000    -0.0013    -0.0010     0.0023     0.0068  0.0006 

10.0000     0.1611     0.1750     0.0077     0.0222  0.0075 

11.0000     0.1031     0.0200     0.0019     0.0055  0.0009 

12.0000    -0.0668    -0.0750     0.0019     0.0056  0.0006 

MMLE    STABILITY  &  CONTROL  DERIVATIVES 

CYJb        Cl_b        CNJb        Cljp        CN_p  Cl_r 

-0.5630    -0.0746     0.0715      -0.4030    -0.0602  0.1029 

CN_r       Clda       CNda       CYdr       Cldr  CNdr 

-0.1189     0.1307    -0.0013       0.1611     0.1031  -0.0668 

INITIAL  INPUT  DERIVATIVES 

-0.6000    -0.1500     0.2000      -0.3500    -0.0500  0.1500 

-0.2000     0.0500    -0.0010       0.1750     0.0200  -0.0750 

"TRUTH  DERIVATIVES"  USED  TO  GENERATE  DATA 

-0.5640    -0.0740     0.0710      -0.4100    -0.0575  0.1070 

-0.1250     0.1340    -0.0035       0.1570     0.1070  -0.0720 
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Figure  B-2  8  UAV  Aileron  and  Rudder  Inputs 
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Figure  B-30  UAV  Yaw  Rate  Response 
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Figure  B-31  UAV  Bank  Angle  Response 
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Figure  B-32  UAV  Roll  Rate  Response 
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UAV  SIMULATED  LATERAL-DIRECTIONAL  RESULTS 
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B.   ACTUAL  TEST  FLIGHT  OUTPUT 
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Figure  B-36  F-14A  Pitch  Rate  Response 
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F-14A  LONGITUDINAL  FLIGHT  TEST  RESULTS 
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b.       T-37 
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Figure  B-39  T-37  Elevator  Input 
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Figure  B-41  T-37  Pitch  Rate  Response 
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T-37     ESTIMATED  AND   MEASURED  G  RESPONSE 
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Figure  B-43  T-37  Normal  Acceleration  Response 


T-37  LONGITUDINAL  FLIGHT  TEST  RESULTS 
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